2. Variable-property corrections

2.1. Inner layer

From several studies in the past decades (Morkovin Reference Morkovin1962; Coleman et al. Reference Coleman, Kim and Moser1995; Huang et al. Reference Huang, Coleman and Bradshaw1995; Patel et al. Reference Patel, Peeters, Boersma and Pecnik2015; Trettel & Larsson Reference Trettel and Larsson2016; Modesti & Pirozzoli Reference Modesti and Pirozzoli2016; Patel, Boersma & Pecnik Reference Patel, Boersma and Pecnik2017; Zhang et al. Reference Zhang, Duan and Choudhari2018), it has been shown that turbulence quantities, when semi-locally scaled (i.e. using $u_\tau ^*$ and $\delta _v^*$ as the relevant velocity and length scales, respectively), collapse well onto their respective incompressible counterparts when plotted as a function of the semi-local coordinate $y^*$ . Some of these quantities, in their classically and semi-locally scaled form are listed in table 1.

Table 1. An example of quantities that are classically and semi-locally scaled in the inner layer. Note that $\bar u^*$ in ${\rm d} \bar u^*/{\rm d}y^*$ represents the semi-locally transformed mean velocity, as defined in (1.3).

Similarly, we argue that if the individual variables collapse when semi-locally scaled, then their model equations written in the semi-locally scaled form must be analogous to those written in the classically scaled form for incompressible flows. Here, we enforce a strict analogy by replacing the individual variables in a classically scaled equation by their semi-locally scaled counterparts.

Let us start by writing the modelled turbulence kinetic energy equation in the inner layer of a canonical incompressible constant-property flow. Neglecting advection terms, the equation reduces to a simple balance between production, dissipation and diffusion of turbulent kinetic energy as

(2.1) \begin{equation} \mu _t \left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 - \beta ^\star \rho _w k \omega + \frac {\rm d}{{\rm d}y}\left [\left (\mu _w + \sigma _k \mu _t\right ) \frac {{\rm d} k}{{\rm d}y} \right ] = 0, \end{equation}

where $\mu _t$ is the eddy viscosity, $\bar u$ the mean velocity, $k$ the turbulence kinetic energy, $\omega$ the specific dissipation rate, $\rho _w$ , $\mu _w$ the density and viscosity at the wall, respectively, and $\beta ^\star$ , $\sigma _k$ are the SST model constants.

Rewriting this equation using the non-dimensional form of the variables, as given in table 1, leads to

(2.2) \begin{equation} {\mu _t^+} \left (\frac {{\rm d}\bar u^+}{{\rm d}y^+}\right )^2 - \beta ^\star k^+ \omega ^+ + \frac {\rm d}{{\rm d}y^+}\left [\left (1 + \sigma _k \mu _t^+\right ) \frac {{\rm d} k^+}{{\rm d}y^+} \right ] = 0, \end{equation}

where the superscript ‘ $+$ ’ denotes the classical wall-based scaling. Now we replace all classically scaled variables with their semi-locally scaled counterparts (refer to table 1), which gives

(2.3) \begin{equation} {\mu _t^*} \left (\frac {{\rm d}\bar u^*}{{\rm d}y^*}\right )^2 - \beta ^\star k^* \omega ^* + \frac {\rm d}{{\rm d}y^*}\left [\left (1 + \sigma _k \mu _t^*\right ) \frac {{\rm d} k^*}{{\rm d}y^*} \right ] = 0, \end{equation}

where the superscript ‘ $*$ ’ denotes semi-local scaling. Rewriting (2.3) using the dimensional form of the variables (see table 1), we get

(2.4) \begin{equation} \frac {\mu _t}{\bar \mu } \left (\frac {\delta _v^*}{u_\tau ^*}\right )^2\left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 - \beta ^\star \frac {k}{{u_\tau ^*}^2} \frac {\omega }{u_\tau ^*/\delta _v^*} + \frac {\rm d}{{\rm d} (y/\delta _v^*)}\left [\left (1 + \sigma _k \frac {\mu _t}{\bar \mu }\right ) \frac {{\rm d }(k/{u_\tau ^*}^2)}{{\rm d} (y/\delta _v^*)} \right ] = 0. \end{equation}

Using the definitions of $u_\tau ^*$ and $\delta _v^*$ from (1.1), we get

(2.5) \begin{equation} \frac {\mu _t}{\bar \mu } \frac {\bar \mu ^2}{\tau _w^2}\left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 - \beta ^\star \frac {\bar \rho }{\tau _w} k \frac {\bar \mu }{\tau _w}\omega + \frac {\rm d}{{\rm d} (y\sqrt {\tau _w\bar \rho }/\bar \mu )}\left [\left (1 + \sigma _k \frac {\mu _t}{\bar \mu }\right ) \frac {{\rm d} (\bar \rho k/\tau _w)}{{\rm d} (y\sqrt {\tau _w\bar \rho }/\bar \mu )} \right ] = 0. \end{equation}

Dividing the equation by $\bar \mu /\tau _w^2$ gives

(2.6) \begin{equation} {\mu _t}\left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 - \beta ^\star {\bar \rho }k \omega + \frac {1}{\bar \mu }\frac {\rm d}{{\rm d} (y\sqrt {\bar \rho }/\bar \mu )}\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {1}{\bar \mu } \frac {{\rm d} (\bar \rho k)}{{\rm d} (y\sqrt {\bar \rho }/\bar \mu )} \right ] = 0. \end{equation}

Applying the chain rule to the diffusion term, we get

(2.7) \begin{equation} {\mu _t}\left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 - \beta ^\star {\bar \rho }k \omega + \frac {1}{\bar \mu }\frac {{\rm d}y}{{\rm d} (y\sqrt {\bar \rho }/\bar \mu )}\frac {{\rm d}}{{\rm d}y}\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {1}{\bar \mu } \frac {{\rm d}y}{{\rm d} (y\sqrt {\bar \rho }/\bar \mu )}\frac {{\rm d} (\bar \rho k)}{{\rm d}y} \right ] = 0, \end{equation}

where ${\rm d} (y\sqrt {\bar \rho }/\bar \mu )/{\rm d}y$ represents the stretching of the wall-normal coordinate due to variable density and viscosity. Note that $y$ in $y\sqrt {\bar \rho }/\bar \mu$ corresponds to the distance from the closest wall. However, $y$ represents this distance only when the origin of the $y$ -axis is located on the wall. To make the equations invariant to the choice of the origin of the $y$ -axis, i.e. to ensure Galilean invariance, we propose to replace $y$ with $\ell$ , where $\ell$ represents the distance from the closest wall.

Finally, introducing a stretching variable $S_y$ , defined as

(2.8) \begin{equation} S_y = \left (\frac {{\rm d} \left (\ell \sqrt {\bar \rho }/\bar \mu \right )}{{\rm d} y}\right )^{-1}=\left (\frac {\sqrt {\bar \rho }}{\bar \mu } + \ell \,\frac {{\rm d} \left (\sqrt {\bar \rho }/\bar \mu \right )}{{\rm d} y}\right )^{-1}, \end{equation}

we get

(2.9) \begin{equation} {\mu _t}\left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 - \beta ^\star {\bar \rho }k \omega + \frac {S_y}{\bar \mu }\frac {\rm d}{{\rm d}y}\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {S_y}{\bar \mu }\frac {{\rm d} (\bar \rho k)}{{\rm d}y} \right ] = 0. \end{equation}

Comparing (2.9) with a standard TKE equation for compressible flows, one can note that the present compressibility corrections only modify the modelling of the diffusion term, just as in the CA/OPDP corrections. It is important to note that since we enforce a strict analogy, these corrections modify the modelling of the total diffusion term, not just the turbulent diffusion term.

For the ease of implementation in existing computational fluid dynamic solvers, these compressibility modifications can be reformulated in the form of a source term ( $\varPhi _k^{\textit{in}}$ ; where the superscript ‘in’ represents inner layer) as follows:

(2.10) \begin{equation} {\mu _t}\left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 - \beta ^\star {\bar \rho }k \omega + \frac {\rm d}{{\rm d} y}\left [\left (\bar \mu + \sigma _k{\mu _t}\right )\frac {{\rm d} k}{{\rm d} y} \right ] + \varPhi ^{\textit{in}}_k= 0, \end{equation}

where the other terms except $\varPhi _k^{\textit{in}}$ represent the standard TKE equation in the inner layer. Then, $\varPhi _k^{\textit{in}}$ can be written as

(2.11) \begin{equation} \varPhi _k^{\textit{in}} = \frac {S_y}{\bar \mu }\frac {\rm d}{{\rm d} y }\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {S_y}{\bar \mu } \frac {{\rm d} (\bar \rho k)}{{\rm d}y} \right ] - \frac {\rm d}{{\rm d} y}\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {{\rm d} k}{{\rm d} y} \right ]. \end{equation}

Here, $\varPhi _k^{\textit{in}} = 0$ for cases with constant mean properties, while $\varPhi _k^{\textit{in}}$ is consistent with the CA/OPDP corrections for cases where $\delta _v^*$ remains uniform in the inner layer.

Repeating the same procedure for the $\omega$ model equation, we obtain

(2.12) \begin{equation} \varPhi _\omega ^{\textit{in}} = \frac {\bar \rho S_y}{\bar \mu ^2}\frac {\rm d}{{\rm d} y }\left [\left (\bar \mu + \sigma _\omega {\mu _t}\right )\frac {S_y}{\bar \mu } \frac {{\rm d} (\bar \mu \omega )}{{\rm d} y} \right ] -\frac {\rm d}{{\rm d} y}\left [\left (\bar \mu + \sigma _\omega \mu _t\right ) \frac {{\rm d} \omega }{{\rm d} y} \right ]. \end{equation}

The standard $k$ - $\omega$ SST model equations along with the source terms defined in (2.11) and (2.12) represent the proposed compressibility corrections in the inner layer.

Figure 1 shows $\mu _t/\bar \mu$ obtained with the $k$ - $\omega$ SST model with (a) no corrections, with (b) CA/OPDP corrections and with (c) the present corrections for numerous compressible turbulent boundary layers (represented by grey lines) from the literature (see § 6 for details on the implementation and § 7, table 2 for details on the cases). The black lines in the figure show an incompressible case of Sillero, Jiménez & Moser (Reference Sillero, Jiménez and Moser2013) at ${Re}_\tau = 1437$ . The collapse on to a single curve in figure 1(c) clearly shows that the present corrections are consistent with the semi-local scaling framework.

Figure 1. Wall-normal distributions of $\mu _t/\bar \mu$ computed using the $k$ - $\omega$ SST model with (a) no corrections, (b) CA/OPDP corrections and (c) present corrections for the zero-pressure-gradient turbulent boundary layers described in § 7 (table 2). The black lines represent the constant-property case of Sillero et al. (Reference Sillero, Jiménez and Moser2013) at ${Re}_\tau =1437$ . For details about the implementation, refer to § 6.

2.2. Outer layer

In the outer layer, the appropriate length scale is $\delta$ (boundary layer thickness; equivalent to $h$ for channels) instead of $\delta _v^*$ . Thus, the corrections discussed so far, which are based on $\delta _v^*$ , are not valid there. To obtain valid compressibility corrections in the outer layer, one may follow the same approach as in the previous section, but using $u_\tau ^*$ and $\delta$ as the relevant scales (Smits & Dussauge Reference Smits and Dussauge2006; Hasan et al. Reference Hasan, Larsson, Pirozzoli and Pecnik2024b ). This leads to (see Appendix A for a detailed derivation)

(2.13) \begin{equation} \varPhi _k^{\textit{out}} = \frac {1}{\sqrt {\bar \rho }}\frac {\partial }{\partial y}\left [\left (\bar \mu +\sigma _{k} \mu _t \right )\frac {1}{\sqrt {\bar \rho }} \frac {\partial (\bar \rho k)}{\partial y} \right ] - \frac {\partial }{\partial y}\left [\left (\bar \mu +\sigma _{k} \mu _t\right ) \frac {\partial k}{\partial y}\right ] \end{equation}

and

(2.14) \begin{equation} \varPhi _\omega ^{\textit{out}} = \frac {\partial }{\partial y}\left [\left (\bar \mu +\sigma _{\omega } \mu _t \right )\frac {1}{\sqrt {\bar \rho }} \frac {\partial (\sqrt {\bar \rho } \omega )}{\partial y} \right ] - \frac {\partial }{\partial y}\left [\left (\bar \mu +\sigma _{\omega } \mu _t\right ) \frac {\partial \omega }{\partial y}\right ]. \end{equation}

Note that these equations are identical to the corrections proposed by Otero Rodriguez et al. (Reference Otero Rodriguez, Patel, Diez Sanhueza and Pecnik2018).

In the SST model, the transport equation for $\omega$ includes a cross-diffusion term, given by

(2.15) \begin{equation} C D_{k \omega }= 2 (1-F_1) \frac {\bar \rho \sigma _{\omega 2}}{\omega } \frac {\partial k}{\partial y} \frac {\partial \omega }{\partial y}, \end{equation}

where $F_1$ is a blending function introduced by Menter (Reference Menter1993) to smoothly transition between the $k$ - $\omega$ and $k$ - $\epsilon$ models. This term is primarily active in the outer region of a boundary layer, where $F_1 \to 0$ . We argue that because the cross-diffusion term originates from the turbulent diffusion term – and given that our proposed corrections target the diffusion term – it is necessary to modify this term for consistency.

Following the same approach as described earlier, the variable-property-corrected cross-diffusion term is given by

(2.16) \begin{equation} C D_{k \omega }= 2 (1-F_1) \frac {\sigma _{\omega 2}}{\sqrt {\bar \rho } \omega } \frac {\partial \bar \rho k}{\partial y} \frac {\partial \sqrt {\bar \rho }\omega }{\partial y}. \end{equation}

Finally, as previously done, we express the corrections on the $CD$ term in the form of a source term as

(2.17) \begin{equation} \varPhi _{CD} = \underbrace {2 (1-F_1) \frac {\sigma _{\omega 2}}{\sqrt {\bar \rho } \omega } \frac {\partial \bar \rho k}{\partial y} \frac {\partial \sqrt {\bar \rho }\omega }{\partial y}}_{\textrm {variable-property-corrected}} - \underbrace {2 (1-F_1) \frac {\bar \rho \sigma _{\omega 2}}{\omega } \frac {\partial k}{\partial y} \frac {\partial \omega }{\partial y}}_{\textrm {conventional}}. \end{equation}

3. Proposed variable-property corrections for the entire boundary layer

Implementing different corrections for the inner and outer layers requires switching between the two forms or blending them with suitable functions. This is essential to preserve the distinct scaling characteristics of these two regions, as demonstrated by Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2024b ) using mixing-length models.

To obtain corrections valid in the entire boundary layer, we propose to blend the inner- ((2.11) and (2.12)) and outer-layer ((2.13) and (2.14)) source terms using a suitable blending function $\mathcal{F}$ as

(3.1) \begin{align} \varPhi _k & = \mathcal{F} \varPhi _k^{\textit{in}} + \left (1-\mathcal{F}\right )\varPhi _k^{\textit{out}} , \end{align}

(3.2) \begin{align} \varPhi _\omega &= \mathcal{F}\varPhi _\omega ^{\textit{in}} + \left (1-\mathcal{F}\right )\varPhi _\omega ^{\textit{out}}, \end{align}

where $\mathcal{F}$ transitions from unity in the viscous sublayer to zero in the outer layer. In this paper, we use $F_1$ as the appropriate blending function; however, other options could be more suitable, whose discussion is deferred to a future study.

With these corrections, the $k$ and $\omega$ transport equations in an SST model are modified as

(3.3) \begin{align} \bar \rho \frac {D k}{Dt} & = \textrm {Conventional TKE terms}+ \varPhi _{k}, \end{align}

(3.4) \begin{align} \bar \rho \frac {D \omega }{Dt} & = \textrm {Conventional} \,\omega \, \textrm {terms} + \varPhi _{\omega }+ \varPhi _{CD}. \end{align}

Lastly, three important points are worth noting. First, only the wall-normal component of the diffusion term in (2.11)–(2.14) is modified; the wall-parallel components remain unchanged. This is because semi-local scaling, and thus the corrections derived from it, are concerned with variations only along the wall-normal direction. However, applying these corrections also in the wall-parallel directions should not influence the results, since mean property gradients along those directions are comparatively negligible. The formulation of the corrections in cases where the wall-normal direction does not align with the coordinate axis (e.g. flow over an inclined wall) is discussed in Appendix B. Second, as expected, the proposed corrections vanish in the free stream region, where the mean properties remain constant. Third, in free-shear flows, these corrections – especially the outer-layer ones – would have minimal impact, leading to little improvement over an uncorrected model (Aupoix Reference Aupoix2004). The inner-layer corrections are irrelevant in such flows due to negligible viscous effects.

4. Intrinsic compressibility corrections

The semi-local scaling approach presented in the previous section accounts for effects associated with mean density and viscosity variations. However, at higher Mach numbers, in addition to variable-property effects, intrinsic compressibility effects also play an important role. These compressibility effects modify the near-wall damping of turbulence, particularly in the buffer layer, resulting in an outward shift in the eddy-viscosity profile and an upward shift in the logarithmic portion of the semi-locally (or TL) transformed mean velocity profile (Hasan et al. Reference Hasan, Larsson, Pirozzoli and Pecnik2023, Reference Hasan, Costa, Larsson, Pirozzoli and Pecnik2025).

To account for this outward shift in the eddy viscosity, Hasan et al. (Reference Hasan, Larsson, Pirozzoli and Pecnik2023) modified the Van Driest damping function in a mixing length model (1.5) as

(4.1) \begin{equation} \mu _t = {\bar \rho u_\tau ^*\kappa yD(y^*,0)} \underbrace {\frac {D(y^*,M_\tau )}{D(y^*,0)}}_{D^{\textit{ic}}}, \end{equation}

where $\bar \rho u_\tau ^* \kappa y D(y^*,0)$ is the semi-local eddy viscosity that accounts for mean property variations alone, and $D^{\textit{ic}}$ is the change in damping caused due to intrinsic compressibility effects. Note that $D^{\textit{ic}}$ is mainly active in the buffer layer and approaches unity in the log-layer and beyond.

Similarly, we propose for the SST model, to multiply the eddy viscosity with a damping function that captures the outward shift due to intrinsic compressibility effects as

(4.2) \begin{equation} \mu _t = {\underbrace {{\frac {a_1 \bar \rho k}{\textit{max}\,(a_1\omega ,F_2\varOmega )}}}_{\textrm {conventional}\,\,\mu _{t}} } ({D^{\textit{ic}}})_{\textit{sst}}, \end{equation}

where the constant $a_1=0.31$ and $\varOmega = {\rm d}\bar u/{\rm d}y$ for canonical flows (Menter Reference Menter1993). The damping function in this equation is defined as

(4.3) \begin{equation} (D^{\textit{ic}})_{\textit{sst}} = \frac {D(R_t,M_t)}{D(R_t, 0)}, \end{equation}

with

(4.4) \begin{equation} D(R_t,M_t) = \left [1 - \mathrm{exp}\left (\frac {-R_t}{K + f(M_t)}\right )\right ]^2, \end{equation}

where $M_t = \sqrt { 2 k}/\bar a$ ( $\bar a$ being the local speed of sound) is the turbulence Mach number and $R_t = \bar \rho k/(\bar \mu \omega )$ is the turbulence Reynolds number. The constant $K$ controls the region in which the damping function $(D^{\textit{ic}})_{\textit{sst}}$ is active (analogous to $A^+$ in (1.5)). To apply the modifications mainly in the buffer layer, we choose $K = 3.5$ .

The function $f(M_t)$ controls how the eddy-viscosity decreases (or shifts outward) with increasing Mach number, analogous to the function $f(M_\tau )$ in (1.5). Based on several high-Mach-number DNS from literature, we find that for the $k$ - $\omega$ SST model,

(4.5) \begin{equation} f(M_t) = 0.39 M_t^ {0.77} \end{equation}

produces accurate results. Refer to Appendix C for details on the formulation and tuning of this function.

5. Modelling the energy equation

We now turn our attention to accurately modelling the energy equation, which provides the density, viscosity and other thermo-physical properties essential for solving the turbulence model equations discussed in the previous sections. We mainly focus on the inner layer, as in this region, the various terms in the energy equation are large in magnitude, which are likely responsible for the high errors observed in the temperature predictions.

The total energy equation in the inner layer of canonical wall-bounded flows can be written as

(5.1) \begin{equation} \frac {\rm d}{{\rm d}y} \left (-\bar {q}_y - \bar {\rho} \! \widetilde{\ v^{\prime \prime }h^{\prime \prime }}\right )+\frac {\rm d}{{\rm d}y} \left (\tilde {u} \bar {\tau }_{xy} -\tilde {u} \bar {\rho} \! \widetilde {\ u^{\prime \prime } v^{\prime \prime }}\right ) +\frac {\rm d}{{\rm d}y}\left (\overline {u_i^{\prime \prime } \tau ^{\prime }_{\textit{ij}}} - \frac {1}{2} \bar {\rho} \!\! \widetilde{\ \ v^{\prime \prime } u_i^{\prime \prime } u_i^{\prime \prime }\ \ \!\!\!}\right ) {- \frac {{\rm d} \bar p}{{\rm d}x}} \tilde {u} = 0, \end{equation}

where $q_y$ is the molecular heat flux in the wall-normal direction, $v$ the wall normal velocity, $h$ the enthalpy and $\tau _{\textit{ij}}$ the viscous shear stress. Here, the tilde denotes density-weighted (Favre) averaging. The first term on the left-hand side represents molecular and turbulent diffusion of enthalpy, the second term represents molecular and turbulent diffusion of mean kinetic energy, and the third term represents molecular and turbulent diffusion of turbulent kinetic energy. The last term represents the work done by the mean pressure gradient in the streamwise direction. For fully developed channel flows, ${\rm d} \bar p/{\rm d}x =-{\tau _w}/{h}$ , and for zero-pressure-gradient boundary layers, ${\rm d} \bar p/{\rm d}x=0$ .

We can simplify the second term on the left-hand side as

(5.2) \begin{equation} \frac {\rm d}{{\rm d}y} \left (\tilde {u} \bar {\tau }_{xy} -\tilde {u} \bar {\rho} \! \widetilde {\ u^{\prime \prime } v^{\prime \prime }}\right ) = \tau _{tot}\frac {{\rm d} \tilde {u}}{{\rm d}y} + \tilde {u}\frac {{\rm d} \tau _{tot}}{{\rm d}y}, \end{equation}

where $\tau _{tot} = \bar {\tau }_{xy} - \bar {\rho } \widetilde {\ u^{\prime \prime } v^{\prime \prime }}$ . For fully developed channel flows, we can write ${{\rm d} \tau _{tot}}/{{\rm d}y} = -{\tau _w}/{h} = {{\rm d}\bar p/{\rm d}x}$ and for ZPG boundary layers, ${{\rm d} \tau _{tot}}/{{\rm d}y} = 0$ in the inner layer. Taking these simplifications into account, we get

(5.3) \begin{equation} \frac {\rm d}{{\rm d}y} \left (-\bar {q}_y - \bar {\rho} \! \widetilde {\ v^{\prime \prime }h^{\prime \prime }}\right )+ (\bar {\tau }_{xy} - \bar {\rho } \!\widetilde {\ u^{\prime \prime } v^{\prime \prime }}) \frac {{\rm d} \tilde {u}}{{\rm d}y} +\frac {\rm d}{{\rm d}y}\left (\overline {u_i^{\prime \prime } \tau ^{\prime }_{\textit{ij}}} - \frac {1}{2} \bar {\rho } \! \widetilde {\ v^{\prime \prime } u_i^{\prime \prime } u_i^{\prime \prime }}\right ) = 0. \end{equation}

The terms in (5.3) are modelled as follows (Wilcox et al. Reference Wilcox2006):

(5.4) \begin{align}& -\frac {{\rm d}\bar {q}_y}{{\rm d}y} = \frac {\rm d}{{\rm d}y}\left (\frac {\bar \mu c_p}{Pr} \frac {{\rm d} \bar T}{{\rm d}y}\right ), \quad -\frac {{\rm d}\bar {\rho } \!\widetilde {\ v^{\prime \prime }h^{\prime \prime }}}{{\rm d}y} = \frac {\rm d}{{\rm d}y}\left (\frac { \mu _t c_p}{Pr_t} \frac {{\rm d} \bar T}{{\rm d}y}\right ), \end{align}

(5.5) \begin{align}&\qquad \tau _{xy}\frac {{\rm d} \tilde {u}}{{\rm d}y} = \bar \mu \left (\frac {{\rm d} \tilde {u}}{{\rm d}y}\right )^2, \quad -\bar {\rho } \!\widetilde {\ u^{\prime \prime } v^{\prime \prime }}\frac {{\rm d} \tilde {u}}{{\rm d}y} = \mu _t \left (\frac {{\rm d} \tilde {u}}{{\rm d}y}\right )^2, \end{align}

where $\textit{Pr}_t$ is the turbulent Prandtl number. The last term on the left-hand side in (5.3) represents the total diffusion of TKE, which is modelled as (see § 2.1, (2.9))

(5.6) \begin{equation} \frac {\rm d}{{\rm d} y}\left (\overline {u_i^{\prime \prime } \tau ^{\prime }_{\textit{ij}}} - \frac {1}{2} \bar {\rho} \! \widetilde {\ v^{\prime \prime } u_i^{\prime \prime } u_i^{\prime \prime }}\right )= \frac {S_y}{\bar \mu }\frac {\rm d}{{\rm d} y }\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {S_y}{\bar \mu } \frac {{\rm d} (\bar \rho k)}{{\rm d}y} \right ]. \end{equation}

With these modelled terms, (5.3) can be written as

(5.7) \begin{equation} \frac {\rm d}{{\rm d}y}\left (\left [\frac {\bar \mu c_p}{Pr}+\frac {\mu _t c_p}{Pr_t}\right ] \frac {{\rm d} \bar T}{{\rm d}y}\right ) = -\left (\bar \mu + \mu _t \right ) \left (\frac {{\rm d} \tilde {u}}{{\rm d}y}\right )^2 - \frac {S_y}{\bar \mu }\frac {\rm d}{{\rm d} y }\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {S_y}{\bar \mu } \frac {{\rm d} (\bar \rho k)}{{\rm d}y} \right ]. \end{equation}

In the inner layer, the total diffusion of TKE (last term in (5.7)) is balanced by its production and dissipation (see (2.9)), such that

(5.8) \begin{equation} \frac {S_y}{\bar \mu }\frac {\rm d}{{\rm d}y}\left [\left (\bar \mu + \sigma _k {\mu _t}\right )\frac {S_y}{\bar \mu } \frac {{\rm d} (\bar \rho k)}{{\rm d}y} \right ] = -{\mu _t}\left (\frac {{\rm d}\bar u}{{\rm d}y}\right )^2 + {\bar \rho } \epsilon . \end{equation}

Using this relation to substitute the total diffusion of TKE in (5.7), we get

(5.9) \begin{equation} \frac {\rm d}{{\rm d}y}\left (\left [\frac {\bar \mu c_p}{Pr}+\frac {\mu _t c_p}{Pr_t}\right ] \frac {{\rm d} \bar T}{{\rm d}y}\right ) = -\left (\bar \mu + \mu _t \right ) \left (\frac {{\rm d} \tilde {u}}{{\rm d}y}\right )^2 + \mu _t \left (\frac {{\rm d} \tilde {u}}{{\rm d}y}\right )^2 - \bar \rho \epsilon , \end{equation}

which further simplifies to

(5.10) \begin{equation} \frac {\rm d}{{\rm d}y}\left (\left [\frac {\bar \mu c_p}{Pr}+\frac {\mu _t c_p}{Pr_t}\right ] \frac {{\rm d} \bar T}{{\rm d}y}\right )= - \underbrace {\bar \mu \left (\frac {{\rm d} \tilde {u}}{{\rm d}y}\right )^2 - \bar \rho \epsilon }_{\varPhi _e}, \end{equation}

where $\varPhi _e$ represents the source term in this equation.

Commonly, the TKE dissipation term in (5.10) is assumed to be equal to the production term in the inner layer, thereby implying equilibrium, i.e. $\bar \rho \epsilon \approx \mu _t ({\rm d}\bar u/{\rm d}y)^2$ (Larsson et al. Reference Larsson, Kawai, Bodart and Bermejo-Moreno2016; Bose & Park Reference Bose and Park2018; Huang et al. Reference Huang, Coleman, Spalart and Yang2023). In other words, the TKE diffusion term in (5.7) is assumed to be zero. However, this assumption breaks down near the wall and is responsible for inaccuracies in predicting the temperature peak in high-speed, cooled-wall boundary layers (see § 7). This implies that accurate modelling of $\epsilon$ (or the TKE diffusion term) in the inner layer is essential for accurate temperature predictions.

Another possible source of error in the energy equation (5.10) is the turbulent Prandtl number, which is often assumed to be a constant and equal to 0.9 (Wilcox et al. Reference Wilcox2006). However, it is well known that for cooled-wall boundary layers, $\textit{Pr}_t$ is not a constant but rather varies substantially in the near-wall region. Furthermore, at the location of the temperature peak, $\textit{Pr}_t$ is undefined (Griffin et al. Reference Griffin, Fu and Moin2023; Chen et al. Reference Chen, Gan and Fu2024). In § 7, we will briefly analyse the sensitivity of our results with respect to different values of $\textit{Pr}_t$ .

6. Implementation

The proposed corrections primarily differ from the CA/OPDP corrections within the inner layer, while they remain identical in the outer layer (except for the cross-diffusion term discussed in § 3). Consequently, for a meaningful comparison between the two approaches, it suffices to focus on solving the inner layer. This also allows us to gauge the effectiveness of the proposed inner-layer corrections without any potential influence from the outer layer.

In this section, we present the implementation of the proposed corrections for simple canonical flows that can be modelled as one-dimensional problems, such as the inner layer of zero-pressure-gradient (ZPG) boundary layers and fully developed channel flows. In these flows, mainly the effect of $\varPhi _k^{\textit{in}}$ , $\varPhi _\omega ^{\textit{in}}$ and the damping function $D^{\textit{ic}}_{\textit{sst}}$ is tested. Subsequently, in § 8, we test the full corrections, i.e. $\varPhi _k^{\textit{in}}$ and $\varPhi _\omega ^{\textit{in}}$ in the inner layer, blended with $\varPhi _k^{\textit{out}}$ and $\varPhi _\omega ^{\textit{out}}$ in the outer layer ((3.1) and (3.2)), together with the damping function $D^{\textit{ic}}_{\textit{sst}}$ and the correction for the cross-diffusion term $\varPhi _{CD}$ .

6.1. Zero-pressure-gradient boundary layers

To solve the inner layer of boundary layers, we follow the methodology outlined in § 4.6.3 of Wilcox et al. (Reference Wilcox2006) for incompressible flows. This approach involves solving the integrated streamwise momentum equation, also known as the total stress balance equation, in conjunction with one-dimensional turbulence model equations to determine the inner-layer velocity profile in a turbulent boundary layer. In this work, we extend this framework to compressible flows.

For compressible boundary layers, the inner-layer velocity and temperature profiles are obtained by solving the integrated forms of the momentum and energy equations, given as

(6.1) \begin{equation} \begin{aligned} \frac {{\rm d} \bar {u}}{{\rm d} y}=&\frac {\tau _w}{\bar {\mu }+{\mu _t}}, \quad \text{and} \\ c_p\left (\frac {\bar \mu }{Pr}+\frac {\mu _t}{Pr_t}\right ) \frac {{\rm d} \bar T}{{\rm d}y} =& -\frac {\mu _w c_p}{Pr}\left (\frac {{\rm d} \bar T}{{\rm d}y}\right )_w - \int _0^{y} \varPhi _e \, {\rm d}y, \end{aligned} \end{equation}

respectively. These equations are solved in a domain spanning from $y=0$ to $y=0.2\delta$ , where $y = 0.2\delta$ is arbitrarily taken to be the edge of the inner layer (Örlü et al. Reference Örlü, Fransson and Alfredsson2010).

The dynamic viscosity $\bar \mu$ in (6.1) is computed using Sutherland’s law, the mean density is obtained as $\bar \rho / \rho _w = T_w/\bar T$ and the Prandtl number $\textit{Pr}$ is equal to 0.72. The first term on the right-hand side of the energy equation, scaled by wall-based quantities ( $\rho _w u_\tau c_p T_w$ ), corresponds to the non-dimensional wall heat flux $B_q$ .

The eddy viscosity $\mu _t$ in (6.1) is computed by solving the one-dimensional SST model (without advection terms), analogous to the one-dimensional $k$ - $\omega$ model described by Wilcox et al. (Reference Wilcox2006) (see their (4.184)). This model is solved with the source terms $\varPhi _k^{\textit{in}}$ and $\varPhi _\omega ^{\textit{in}}$ described in (2.11) and (2.12), along with the standard eddy viscosity formulation being multiplied by $(D^{\textit{ic}})_{\textit{sst}}$ (4.3). The boundary conditions for the SST model are defined as

(6.2) \begin{equation} \begin{array}{ccccc} k = 0, & \omega = \dfrac {60 \bar \mu }{\beta _1 \bar \rho {(\Delta y)}^2} & \text{ at } & y = 0, & \text{and}\\ k = \dfrac {{u_\tau ^*}^2}{\sqrt {\beta ^\star }}, & \omega = \dfrac {{u_\tau ^*}}{\sqrt {\beta ^\star } \kappa y} & \text{ at } & y = 0.2\delta , & \end{array} \end{equation}

where $\Delta y$ is the distance to the next grid point away from the wall. The first row of (6.2) corresponds to the wall boundary condition described by Menter (Reference Menter1993), and the second row corresponds to the log-layer asymptotic solution for $k$ and $\omega$ as described by Wilcox et al. (Reference Wilcox2006) (see their (4.185)), but adapted to compressible flows.

The turbulent Prandtl number $\textit{Pr}_t$ in (6.1) is assumed to be a constant and equal to 0.9. However, to show the sensitivity of the results with respect to the constant $\textit{Pr}_t$ assumption, we will also compute some of the results with $\textit{Pr}_t$ from DNS and then compare them with $\textit{Pr}_t=0.9$ .

Lastly, the source term $\varPhi _e$ in (6.1) was defined in (5.10) as

(6.3) \begin{equation} \varPhi _e = \bar \mu ({\rm d} \bar u/{\rm d}y)^2 + \bar \rho \epsilon . \end{equation}

As discussed earlier, a common way to model the TKE dissipation rate is to assume it to be equal to the production term ( $\epsilon =\mu _t ({\rm d} \bar u/{\rm d}y)^2$ ), such that the source term becomes

(6.4) \begin{equation} \varPhi _{e,1} = \bar \mu ({\rm d} \bar u/{\rm d}y)^2 + \mu _t ({\rm d} \bar u/{\rm d}y)^2, \end{equation}

where ‘ $1$ ’ in the subscript corresponds to the first type of approximation that we will use in this paper. This approximation is inaccurate near the wall, since close to the wall, the production term tends to zero, whereas the dissipation term is finite and balances the total diffusion of TKE.

Therefore, we seek a more accurate representation of $\epsilon$ . One way is to use the TKE dissipation rate estimated by the model itself, which, in the case of SST, is equal to $\epsilon _{ {sst}} = \beta ^\star k \omega$ . However, just like the production term, $\epsilon _{\textit{sst}}$ also approaches zero at the wall, as $k \to 0$ while $\omega$ remains finite (see wall-boundary conditions in 6.2). Consequently, using $\epsilon$ from the SST model would yield results similar to those obtained with $\varPhi _{e,1}$ . To address this inherent limitation, we follow the approach of Rahman, Wallin & Siikonen (Reference Rahman, Wallin and Siikonen2012) and model the dissipation rate as

(6.5) \begin{align} \epsilon _{\textit{eff}} = \sqrt {\epsilon _{ {sst}}^2 + \epsilon _w^2}, \end{align}

where $\epsilon _w$ ensures non-zero dissipation rate at the wall, given by $\epsilon _w = 2 A_\epsilon ({\bar \mu }/{\bar \rho }) ({{\rm d} \bar u}/{{\rm d}y})^2$ , with $A_\epsilon = 0.09$ (Rahman et al. Reference Rahman, Wallin and Siikonen2012), and where $\epsilon _{ {eff}}$ represents the effective dissipation rate. With this, the second approximation of $\varPhi _e$ that we will use in this paper is given as

(6.6) \begin{equation} \varPhi _{e,2} = \bar \mu ({\rm d} \bar u/{\rm d}y)^2 + \bar \rho \epsilon _{\textit{eff}}. \end{equation}

We solve the above-mentioned equations iteratively until convergence; see Hasan et al. (Reference Hasan, Elias, Menter and Pecnik2024a ) for more details on the solver.

6.2. Fully developed channel flows

For channel flows, instead of solving only for the inner layer, we solve for the entire domain (spanning from $y=0$ to $y=2h$ , $h$ being the channel half-height). However, we do not switch to the outer-layer compressibility corrections, but consistently use the inner-layer corrections in the entire domain. This is because the velocity profile in the outer layer of channel flows closely follows the logarithmic profile of the inner layer, meaning the error introduced by not solving the outer layer with the correct corrections is minimal.

For fully developed channel flows, we solve the momentum and energy equations, given as

(6.7) \begin{equation} \begin{aligned} \frac {\rm d}{{\rm d} y}\left [\left ({\bar \mu }+\mu _t\right ) \frac {{\rm d} \bar u}{{\rm d} y}\right ] & = - \frac {\tau _w}{h}, & \quad & \frac {\rm d}{{\rm d}y}\left (\left [\frac {\bar \mu c_p}{Pr}+\frac {\mu _t c_p}{Pr_t}\right ] \frac {{\rm d} \bar T}{{\rm d}y}\right )= -\varPhi _e, \end{aligned} \end{equation}

where $\bar \mu$ is computed using a power-law $(\bar T/T_w)^n$ with an exponent of $0.75$ for high-Mach-number flows and $0.7$ for the low-Mach-number cases. The mean density is computed as $\bar \rho /\rho _w = T_w/\bar T$ , and the Prandtl number $\textit{Pr}$ is assumed to be a constant and equal to $0.72$ for the high-Mach-number cases, and is defined as $Pr = (\bar T/T_w)^{0.7}$ for the low-Mach-number cases. The different choices for $\bar {\mu }$ and $\textit{Pr}$ are consistent with those used in the respective DNSs of these cases.

The eddy-viscosity $\mu _t$ in (6.7) is obtained by solving the one-dimensional SST model with the proposed corrections (similar to ZPG boundary layers discussed previously), along with the wall boundary conditions ( $k=0$ and $\omega = {60 \bar \mu }/[\beta _1 \bar \rho (\Delta y)^2]$ ) imposed on the two walls located at $y=0$ and $y=2h$ .

The turbulent Prandtl number ( $\textit{Pr}_t$ ) in (6.7) is assumed to be constant, with $\textit{Pr}_t = 0.9$ for high-Mach-number channel flows, similar to ZPG boundary layers, and $\textit{Pr}_t = 1$ for low-Mach-number flows. The higher value for the low-Mach-number cases is based on the results of Patel et al. (Reference Patel, Boersma and Pecnik2017), where $\textit{Pr}_t$ is reported to be close to 1, or even slightly higher for these cases. This is likely due to the strong similarity between momentum and energy equations in these flows.

The source term $\varPhi _e$ in (6.7) for high-Mach-number channel flows can be approximated as $\varPhi _{e,1}$ or $\varPhi _{e,2}$ , similar to ZPG boundary layers. In contrast, for low-Mach-number flows, the dissipation of mean and turbulent kinetic energy does not contribute to the right-hand side of (6.7). This can be explained as follows: $\varPhi _{e,1}$ (or $\varPhi _{e,2}$ ) scales with $\rho _{r} \mathcal{U}^3/\mathcal{L}$ , where $\rho _{r}$ is a reference density scale, and $\mathcal U$ and $\mathcal{L}$ are relevant velocity and length scales. respectively. However, the left-hand side of the energy (6.7) scales with $\rho _{r} \mathcal{U} c_{p,r} T_{r}/\mathcal{L}$ , where $c_{p,r}$ and $T_{r}$ are reference scales for $c_p$ and $\bar T$ , respectively. The ratio of these scales corresponds to an Eckert number, defined as $\mathcal{U}^2/(c_{p,r} T_r)$ , which is proportional to the square of the Mach number. For flows at low (or zero) Mach numbers, this ratio tends to zero which explains why $\varPhi _{e,1}$ and $\varPhi _{e,2}$ are ineffective in (6.7). Thus, in those flows, temperature and hence property variations are created by adding a user-defined heat source which is uniform in the domain. For the two low-Mach-number gas-like cases considered here, this external heat source is equal to $17.55 (\mu _w c_p/ Pr_w) T_w/h$ (Patel et al. Reference Patel, Peeters, Boersma and Pecnik2015) and $75 (\mu _w c_p/ Pr_w) T_w/h$ (Pecnik & Patel Reference Pecnik and Patel2017).

All the equations described here are solved in their wall-scaled form. For more details on the solver, refer to the jupyter notebook (Hasan et al. Reference Hasan, Elias, Menter and Pecnik2024a ).

8. Testing the proposed corrections on a two-dimensional boundary layer

In the previous sections, we discussed the implementation and results for the inner layer of ZPG boundary layers and channel flows, which can be modelled using one-dimensional equations. However, it is so far unclear how the proposed corrections perform in a two-dimensional boundary layer simulation, where the outer layer is also solved for.

To address this point, we performed a flat plate boundary layer simulation for the $M_\infty = 10.9$ , $T_w/T_r = 0.2$ case of Huang et al. (Reference Huang, Duan and Choudhari2022), using Ansys Fluent. The computational domain spans 1.5 m in the streamwise direction and 0.2 m in the wall-normal direction, discretised using a structured mesh consisting of 1600 cells along the streamwise direction and 300 cells in the wall-normal direction. The mesh is stretched in the wall-normal direction such that the first grid point lies within one viscous wall unit from the wall. With this mesh, we achieve a resolution of $\Delta x^+ \approx 97$ and $\Delta y^+_{{min}} \approx 0.014$ at the streamwise location whose profiles are shown in figure 10. This resolution is finer than that adopted by Danis & Durbin (Reference Danis and Durbin2022) for the same case.

Figure 10. Computed (a) full mean velocity and (b) temperature profiles compared with the DNS (coloured solid lines) for the following turbulent boundary layer: $M_\infty =10.9$ , $T_w/T_r=0.2$ , ${Re}_\tau =774$ (Huang et al. Reference Huang, Duan and Choudhari2022). Refer to the legend for line types. ‘Baseline’ represents the SST model without corrections, while ‘Present’ signifies the corrections proposed in this paper ((2.17), (3.1), (3.2) and (4.3)). Note that, in contrast with figure 3, where only the inner layer was solved, here, both the baseline and present results are obtained by solving the entire boundary layer. Also note that these results are computed with (8.1) as the energy equation, except that for the baseline case, $\varPhi _k=0$ . The solid black lines correspond to the profiles obtained after solving the inner-layer equations (6.1) with $\varPhi _e = \bar \mu ({\rm d}\bar u/{\rm d}y)^2 + \bar \rho \epsilon _{\textit{sst}}$ .

Figure 11. (a) Skin-friction and (b) heat transfer coefficients compared with the DNS (symbols) for the $M_\infty =10.9$ , $T_w/T_r=0.2$ turbulent boundary layer of Huang et al. (Reference Huang, Duan and Choudhari2022). The dotted and solid grey lines correspond to the quantities obtained with the uncorrected and corrected SST models, respectively. The black vertical lines indicate an error margin of +/-5 %. Note that in panel (b), the $y$ -axis has been scaled by a factor of 1.1, corresponding to the assumed value of the ratio $2\,c_h / c_f$ .

A velocity-inlet boundary condition is applied at the inlet, prescribing a constant streamwise velocity of 1778.4 m s⁻¹, pressure of 1966.1 Pa and temperature of 64.4 K. At the wall, a no-slip boundary condition is imposed along with a constant wall temperature of 300 K. A pressure far-field boundary condition is specified at the top boundary, with the pressure and temperature set to 1966.1 Pa and 64.4 K, matching their respective values at the inlet. At the outlet, a pressure-outlet boundary condition is used, prescribing the pressure to be equal to the free stream pressure, and the backflow temperature is set to 300 K. Note that the free stream temperature specified in the DNS is 66.5 K; here, we slightly reduce this value to account for the temperature jump across the oblique shockwave originating from the leading-edge, such that the post-shock temperature is approximately 66.5 K. (The pressure also changes across the shockwave; however, the ratio of its post-shock value to the wall value remains unchanged. In contrast, because the wall temperature is fixed, the corresponding ratio for temperature is affected. This, in turn, alters the density ratio and ultimately influences the computed skin friction and heat transfer coefficients in figure 11.) With these boundary conditions, we perform the simulations using Fluent’s density based solver. Further details on the simulation set-up (for instance, the case file) can be found on our GitHub repository (Hasan et al. Reference Hasan, Elias, Menter and Pecnik2024a ).

To implement the proposed corrections in Fluent, we use user-defined functions (UDFs). In these UDFs, we compute the source terms $\varPhi _{k}$ , $\varPhi _{\omega }$ and $\varPhi _{CD}$ , and add them on the right-hand side of the conventional SST model as indicated in (3.3) and (3.4). We also redefine the eddy-viscosity as in (4.2). The reader is referred to our GitHub repository (Hasan et al. Reference Hasan, Elias, Menter and Pecnik2024a ) for more details on the UDFs used to generate the results shown in the following.

To obtain the mean temperature profile, we solve the full energy equation:

(8.1) \begin{align} \frac {\partial \left [\bar \rho (e + \widetilde {u}_i\widetilde {u}_i + k)\right ]}{\partial t}+\frac {\partial \left [\bar \rho \widetilde {u}_{\!j} (h + \widetilde {u}_i\widetilde {u}_i + k)\right ]}{\partial x_{\!j}} &= \frac {\partial }{\partial x_{\!j}}\left [\left (\frac {\bar {\mu }c_p}{P r}+\frac {\mu _t c_p}{P r_t}\right ) \frac {\partial \bar {T}}{\partial x_{\!j}}\right ] \nonumber \\ &\quad+\frac {\partial }{\partial x_{\!j}}\left [\widetilde {u}_i\left (2\left (\bar {\mu }+\mu _t\right ) \widetilde {S}_{\textit{ij}}\right )\right ] \nonumber \\ &\quad+ \underbrace {\frac {\partial }{\partial x_{\!j}}\left [\left (\bar \mu +\sigma _{k} \mu _t\right ) \frac {\partial k}{\partial x_{\!j}}\right ] + \varPhi _k}_{\text{variable-property-corrected}\,(\S\, 2.1)}, \end{align}

where $e$ and $h$ are the internal energy and enthalpy per unit mass, and $S_{\textit{ij}}$ is the strain rate tensor. (Refer to § 5 for a discussion on the energy equation in the inner layer.) Adding $\varPhi _k$ to the conventional equation implies that the viscous and turbulent diffusion of TKE are modelled after accounting for variable-property effects, consistent with the diffusion term in the TKE transport equation (see § 2.1).

Note that in Fluent, TKE is not included in the definition of total energy, and its diffusion does not appear in the energy equation. Therefore, to correctly solve (8.1), it is necessary to add the TKE transport term $\bar {\rho }\, {\rm D}k/{\rm D}t$ to the left-hand side of the energy equation and the TKE diffusion term to the right-hand side. We achieve this in Fluent by introducing a source term on the right-hand side of the energy equation, equal to the TKE diffusion term minus $\bar {\rho }\, {\rm D}k/{\rm D}t$ . This term is effectively equal to $-P_k + \bar {\rho }\,\epsilon$ , where $P_k$ denotes the production of TKE. For further details, refer to the UDF provided by Hasan et al. (Reference Hasan, Elias, Menter and Pecnik2024a ).

In the inner layer, (8.1) reduces to (5.7). Interestingly, with $k$ taken from the SST model, solving (8.1) in the inner layer is equivalent to solving (6.1) with $\varPhi _e = \bar \mu (\mathrm{d} \bar u/\mathrm{d}y)^2 + \bar \rho \epsilon _{\textit{sst}}$ . However, as discussed in (6.5), solving such an equation in the inner layer would lead to inaccuracies in predicting the mean temperature, since $\epsilon _{\textit{sst}}$ tends to zero at the wall, instead of a non-zero value. Thus, there, we proposed replacing $\epsilon _{\textit{sst}}$ with $\epsilon _{\textit{eff}}$ (6.5), or in other words, we proposed adding a source term, namely, $\bar \rho \epsilon _{\textit{eff}} - \bar \rho \epsilon _{\textit{sst}}$ to $\bar \mu ({\rm d} \bar u/{\rm d}y)^2 + \bar \rho \epsilon _{\textit{sst}}$ , such that effectively we arrive at $\varPhi _{e,2}$ , i.e. (6.6).

In the same spirit, to include the effects of $\epsilon _w$ in the two-dimensional solver, one needs to add a source term equal to $\bar \rho \epsilon _{\textit{eff}} - \bar \rho \epsilon _{\textit{sst}}$ on the right-hand side of (8.1). However, including such a source term does not lead to any improvement in the temperature profile (not shown), unlike what one would expect based on the one-dimensional results, where including the effects of $\epsilon _w$ improved the temperature profiles (see grey solid lines in figure 3 b). This is explained as follows: in the one-dimensional solver (6.1), the wall heat flux (first term on the right-hand side) is fixed and provided as a boundary condition. Thus, any increase in $\varPhi _e$ resulting from adding the source term $\bar \rho \epsilon _{\textit{eff}} - \bar \rho \epsilon _{\textit{sst}}$ would directly reduce the conductive heat flux, thereby reducing the temperature peak. Contrastingly, in the two-dimensional solver, the wall temperature is fixed, and the heat flux is allowed to vary. Now, since the integral of the added source term from the wall to the free stream at a particular streamwise location is non-zero, it directly results in higher wall heat flux, rather than a lower conductive flux close to the wall. Thus, the temperature profile does not improve. Providing an accurate source term using the SST model, without erroneously adding to the wall heat flux, still remains unclear. Nevertheless, here, we provide the results using (8.1), i.e. without accounting for the influence of $\epsilon _w$ .

Figures 10(a) and 10(b) show the mean velocity and temperature profiles, respectively, for the $M_\infty =10.9$ case described previously. These profiles are extracted from the streamwise location where the friction Reynolds number matches the DNS value, i.e. ${Re}_\tau \approx 774$ . The solid coloured lines represent the DNS profiles, while the dotted grey lines depict the predictions from the conventional SST model. The solid grey lines show the results obtained using the corrected SST model (3.3) and (3.4), along with the modified eddy viscosity (4.2). The solid black lines correspond to the results obtained after solving the one-dimensional equations with the present inner-layer corrections and with $\varPhi _e = \bar \mu ({\rm d} \bar u/{\rm d}y)^2 + \bar \rho \epsilon _{\text{sst}}$ (see § 6.1). This $\varPhi _e$ is consistent with equation (8.1) in the inner layer, as discussed earlier.

As seen in figure 10, including the proposed corrections significantly improves the velocity profile compared with the one obtained with the baseline SST model (compare the dotted and solid grey lines). This improvement mainly originates from the inner-layer corrections. We confirmed this by setting $\varPhi _k^{\textit{out}}$ , $\varPhi _\omega ^{\textit{out}}$ and $\varPhi _{CD}$ equal to zero in our simulations, such that we solve the baseline model in the outer layer. The resulting profiles were similar to those obtained when the outer-layer corrections were included. In contrast to the velocity, the temperature profiles remain largely unchanged.

Figure 10 also shows that, in the inner layer, the velocity and temperature profiles computed using Fluent (grey solid line) collapses well on to the respective profiles computed using the one-dimensional approach discussed earlier (black solid lines). Note that in Fluent, the molecular Prandtl number for air is close to 0.74, and the turbulent Prandtl number is assumed to be 0.85. Thus, to be consistent, the one-dimensional results are also computed with these $\textit{Pr}$ and $\textit{Pr}_t$ values.

Finally, in figure 11, we present the distribution of the skin-friction ( $c_f = \tau _w/[0.5 \rho _\infty U_\infty ^2]$ ) and heat transfer ( $c_h = q_w/[\rho _\infty U_\infty c_p (T_r-T_w)]$ ) coefficients along the plate, as a function of ${Re}_\tau$ . The filled symbols represent the DNS data of Huang et al. (Reference Huang, Duan and Choudhari2022), and the black vertical lines represent an error margin of +/-5 %. Clearly, the $c_f$ and $c_h$ predictions obtained from the corrected SST model are closer to the DNS values than those obtained from the conventional model (compare the dotted and solid grey lines).