Near-wall model for compressible turbulent boundary layers based on an inverse velocity transformation

Abstract

In this work, a near-wall model, which couples the inverse of a recently developed compressible velocity transformation (Griffin et al., Proc. Natl Acad. Sci., vol. 118, 2021, p. 34) and an algebraic temperature–velocity relation, is developed for high-speed turbulent boundary layers. As input, the model requires the mean flow state at one wall-normal height in the inner layer of the boundary layer and at the boundary-layer edge. As output, the model can predict mean temperature and velocity profiles across the entire inner layer, as well as the wall shear stress and heat flux. The model is tested in an a priori sense using a wide database of direct numerical simulation high-Mach-number turbulent channel flows, pipe flows and boundary layers (48 cases, with edge Mach numbers in the range 0.77–11, and semi-local friction Reynolds numbers in the range 170–5700). The present model is significantly more accurate than the classical ordinary differential equation (ODE) model for all cases tested. The model is deployed as a wall model for large-eddy simulations in channel flows with bulk Mach numbers in the range 0.7–4 and friction Reynolds numbers in the range 320–1800. When compared to the classical framework, in the a posteriori sense, the present method greatly improves the predicted heat flux, wall stress, and temperature and velocity profiles, especially in cases with strong heat transfer. In addition, the present model solves one ODE instead of two, and has a computational cost and implementation complexity similar to that of the commonly used ODE model.

1. Introduction

The largest driver of computational cost in numerical simulations of wall-bounded turbulence is typically the numerical resolution in the near-wall region. In scale-resolving simulations, e.g. wall-resolved (WR) large-eddy simulations (LES), high spatial and temporal resolutions are required to simulate accurately the small-scale eddies near walls. Wall models, or approximate boundary conditions, can be employed to reduce the near-wall resolution requirements. The computational cost (the number of grid points multiplied by the number of time steps) for the simulation of a turbulent boundary layer scales with the Reynolds number as $Re^{2.7}$ for WRLES and $Re^{1.1}$ for wall-modelled (WM) LES (Yang & Griffin 2021). Thus wall models lead to substantial cost savings for high-Reynolds-number applications. In simulations of the Reynolds-averaged Navier–Stokes (RANS) equations, high spatial resolution is also required to resolve the steep near-wall gradients in the mean flow. Therefore, wall models – typically referred to as wall functions in the RANS context – can also greatly accelerate numerical simulations.

The present work focuses on the paradigm of wall-stress modelling (Larsson et al. 2016; Bose & Park 2018) for LES. These models were derived from RANS analysis of boundary layers, and typically invoke a zero-equation RANS model such as the Prandtl mixing length argument (Prandtl 1925), which models the turbulence length scale as a linear function of the wall-normal distance. An empirical damping function is introduced following Van Driest (1956) to ensure the correct near-wall scaling of the mixing length. RANS models have naturally been used widely as boundary conditions for under-resolved RANS simulations (e.g. Abrahamson & Brower 1988; Lien, Kalitzin & Durbin 1998; Goncalves & Houdeville 2001; Parente et al. 2011). In this context, such a model is typically referred to as a wall function. Cabot (1995) and Cabot & Moin (2000) showed that the mixing length RANS model is suitable for use as a boundary condition for the LES equations, i.e. for deployment as a wall-stress model. Specifically, they invoke the one-dimensional simplification of the RANS streamwise momentum equation. That is,

(1.1)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}y}\left((\bar{\mu} + \bar{\mu}_t)\, \frac{\mathrm{d}\tilde{U}}{\mathrm{d}y} \right) = 0, \end{equation}

where $\bar {\mu }$, $\bar {\mu }_t$ and $\tilde {U}$ are the molecular dynamic viscosity, eddy viscosity and velocity profiles, respectively, and $y$ is the wall-normal coordinate. Here, $\overline {({\cdot })}$ denotes the Reynolds average, and $\widetilde {({\cdot })}$ denotes the Favre (density-weighted) average. Throughout this work the Favre (density-weighted) averaged RANS and LES equations are employed. The eddy viscosity is further modelled as

(1.2)\begin{equation} \bar{\mu}_t = \kappa y \bar{\rho}\, \sqrt{\tau_w/\bar{\rho}} ( 1 - \exp (y^+{/}A^+) )^2, \end{equation}

where $\bar {\rho }(y)$ is the density profile. The subscript $w$ denotes quantities evaluated at the wall, and $\tau _w = \bar {\mu }_w ({\rm d}\tilde {U}/{{\rm d} y})_w$ is the wall shear stress. The superscript $+$ denotes non-dimensionalization by the friction velocity $u_\tau = \sqrt {\tau _w/\bar {\rho }_w}$, $\bar {\rho }_w$ and the kinematic wall viscosity $\bar {\nu }_w=\bar {\mu }_w/\bar {\rho }_w$.

The von Kármán constant $\kappa = 0.41$ and the eddy-viscosity damping coefficient $A^+ = 17$ are adopted following Cabot & Moin (2000).

For an incompressible flow, the density and molecular dynamic viscosity are known constants. In the context of WMLES, the ordinary differential equation (ODE) in (1.1) is solved with two boundary conditions: (1) the no-slip wall condition, and (2) a velocity sample, which is taken from the LES at a wall-normal distance referred to as the matching location. Note that the solution procedure is iterative because the eddy viscosity depends on the wall stress (see (1.2)). The computed wall stress $\tau _w$ is then applied as a momentum-flux boundary condition for the outer LES solver, which completes the two-way coupling of the wall model (inner) solution and the partial differential equation (PDE) (outer) simulation.

For compressible flow, the RANS equation for temperature can be simplified similarly to the one-dimensional form (Larsson et al. 2016; Bose & Park 2018), which results in a second, coupled ODE for the temperature profile, i.e.

(1.3) \begin{equation} \frac{\mathrm{d}}{\mathrm{d}y}\biggl((\bar{\mu} + \bar{\mu}_t) \tilde{U}\,\frac{\mathrm{d}\tilde{U}}{\mathrm{d}y} + C_p \left(\frac{\bar{\mu}}{Pr}+\frac{\bar{\mu}_t}{Pr_t}\right) \frac{\mathrm{d}\tilde{T}}{\mathrm{d}y}\biggr)= 0, \end{equation}

where $\tilde {T}$ is the temperature profile, $C_p$ is the specific heat capacity at constant pressure, $Pr$ is the Prandtl number, and $Pr_t$ is the turbulent Prandtl number, which is assumed to be 0.9 (Larsson et al. 2016). The dependence of molecular dynamic viscosity on temperature can be assumed to follow a power law or Sutherland's law. The ideal gas equation of state closes the system, and the thin boundary-layer assumption implies that the pressure is constant across the inner layer.

In WMLES, the temperature ODE in (1.3) is solved with two additional boundary conditions: (1) the wall temperature, and (2) the temperature at the matching location. Note that the solution procedure is also iterative in that the temperature depends on the velocity solution. The velocity also depends on the temperature through the density and viscosity. Solving two coupled boundary-value problems iteratively introduces a higher degree of nonlinearity compared to the incompressible case, and can prove difficult to converge in flows with strong temperature gradients (strong heat transfer), e.g. as was reported in Fu et al. (2021). In addition to the numerical difficulties, the accuracy of this wall model degrades substantially in flows with strong heat transfer (as will be demonstrated herein).

Improved results for high-speed wall-bounded turbulent flows over cold walls have been obtained by using the semi-local scaling in the damping function (Yang & Lv 2018; Fu, Bose & Moin 2022); however, Iyer & Malik (2019) report that in adiabatic walls, the classical scaling (consistent with the van Driest transformation) is more accurate. This motivates using a recently developed compressible velocity transformation that is accurate for both diabatic and adiabatic turbulent boundary layers (Griffin, Fu & Moin 2021c).

In this work, a wall model for high-speed wall-bounded turbulent flows is developed in § 2. The model is evaluated via a priori testing in § 3, and via a posteriori validation in § 4. Conclusions are drawn in § 5.

2. Model development

There are two principal differences between the present model and the classical ODE-based wall model ((1.1)–(1.3)). (1) Rather than solving an ODE for the compressible velocity profile directly, the incompressible ODE (with constant density and viscosity) is solved, and an inverse compressibility transformation (Griffin et al. 2021c) is employed. (2) Rather than employing a RANS equation for temperature and assuming a constant $Pr_t$, an algebraic temperature–velocity relation is adopted, thus obviating the need to solve a second ODE.

2.1. Inverse compressible velocity transformation

A compressible velocity transformation seeks to map the local mean strain rate of the variable-property compressible flow, ${\rm d}\tilde {U}/{{\rm d} y}$, to the non-dimensional mean strain rate of a constant-property incompressible flow at an equivalent Reynolds number. Upon integration, the transformation maps the compressible velocity profile to an incompressible velocity profile. In this way, a successful transformation can collapse profiles with different Mach numbers and thermal boundary conditions to a single incompressible law of the wall. Coupled with the incompressible profile implied by (1.1), an inverse velocity transformation can recover the compressible velocity profile.

The total-stress-based compressible velocity transformation of Griffin et al. (2021c) is used in this work since it is shown to be accurate in a wide range of flows, including boundary layers with strong heat transfer. This transformation uses the viscous scaling arguments of Trettel & Larsson (2016) and Patel, Boersma & Pecnik (2016) in the near-wall viscous region, and uses a modified version of the turbulence equilibrium arguments of Zhang et al. (2012) for the logarithmic region. The transformation is an algebraic function that relates the local mean strain rate of the compressible flow, ${\rm d}\tilde {U}/{{\rm d} y}$, to the non-dimensional incompressible mean strain rate $S_t^+$ at the same semi-local friction Reynolds number $Re_\tau ^*$, according to the relation

(2.1)\begin{equation} S_t^+{=} \frac{S_{eq}^+}{1+S_{eq}^+{-}S_{TL}^+}, \end{equation}

where $S_{eq}^+=(1/\bar {\mu }^+)\,{\rm d}\tilde {U}^+/{{\rm d} y}^*$ and $S_{TL}^+=\bar {\mu }^+\, {\rm d}\tilde {U}^+/{{\rm d} y}^+$. The superscript $*$ denotes non-dimensionalization by the local density $\rho (y)$, local molecular dynamic viscosity $\mu (y)$, and semi-local friction velocity $u_{sl}=\sqrt {\tau _w/\bar {\rho }(y)}$ (Coleman, Kim & Moser 1995; Huang, Coleman & Bradshaw 1995). The semi-local friction Reynolds number is thus defined as $Re_\tau ^* = \bar {\rho }_e u_{sl} \delta / \bar {\mu }_e$, where the subscript $e$ denotes quantities evaluated at the boundary-layer edge. (Throughout this work, $\delta$ denotes the channel half-height or the boundary-layer thickness.) Note that all variables of the form $S_{({\cdot })}^+$ represent different local non-dimensionalizations of the compressible strain rate, which were designed in prior works with the target of equalling the strain rate implied by the incompressible law of the wall. For example, although $S_{TL}^+$ is equivalent to the viscous stress, it is also a non-dimensionalization of the mean strain rate in a compressible flow; it will recover exactly the incompressible strain rate of a flow with the equivalent viscous stress as long as the compressible flow also obeys $\mu ^+=1$. Additionally, note that the transformation in (2.1) assumes a constant stress layer in the buffer region of the boundary layer, where there is a transition between the underlying viscous and equilibrium transformations. Griffin et al. (2021c) verifies that the deployment of this assumption does not significantly affect the accuracy of the transformation in equilibrium flows, and Bai, Griffin & Fu (2022) verifies the same for boundary layers with moderate pressure gradients.

The inverse velocity transformation is obtained readily by rearranging the transformation algebraically to find

(2.2)\begin{equation} \frac{\mathrm{d}\tilde{U}^+}{\mathrm{d}y^*} = \left(\frac{1}{\bar{\mu}^+ S^+_{t}} - \frac{1}{\bar{\mu}^+} +\sqrt{\bar{\rho}^+} \left(1 + \frac{1}{2\bar{\rho}^+}\,\frac{\mathrm{d}\bar{\rho}^+}{\mathrm{d}y^+}\, y^+{-} \frac{1}{\bar{\mu}^+}\,\frac{\mathrm{d}\bar{\mu}^+}{\mathrm{d}y^+}\,y^+ \right) \right)^{{-}1}. \end{equation}

The incompressible mean strain rate $S_t^+$ is available algebraically from the constant-property version of (1.1), i.e. $\bar {\rho }=\bar {\rho }_w$ and $\bar {\mu }=\bar {\mu }_w$. The incompressible model constants $\kappa$ and $B$ are determined using the aforementioned calibration, but $Re_\tau ^*$ is used in place of $Re_\tau$ since the former is invariant under the velocity transformation. Integrating (2.2) with variable properties yields the targeted compressible velocity profile; the properties are functions of temperature, which will be discussed next.

2.2. Algebraic temperature–velocity relation

In order to close the velocity equation (2.2), the temperature profile must be determined. The classical model uses the constant turbulent Prandtl number assumption to develop a coupled ODE for temperature (1.3). However, the constant Prandtl number assumption has been shown to be less accurate than invoking the generalized Reynolds analogy (GRA) (Zhang et al. 2014). Thus the presently proposed wall model leverages the GRA instead.

The analogy between the conservation equations for momentum and energy has led to the derivation of several algebraic relations between temperature and velocity. Walz's equation (Walz 1969) (also known as the modified Crocco–Busemann relation Busemann 1931; Crocco 1932) leverages the analogy between the conservation equations for momentum and energy to arrive at an algebraic relation between mean temperature and velocity. This relation accounts for non-unity $Pr$ effects via a recovery factor, which is taken as $r=(Pr)^{1/3}$. While this relation is accurate in high-speed adiabatic boundary layers, Duan & Martín (2011) observed that the accuracy degrades significantly in boundary layers with wall heat transfer, and proposed a semi-empirical correction to the relation. Subsequently, this was recast in terms of a GRA (Zhang et al. 2014), thereby introducing the Reynolds analogy factor $s$, which they chose as $s = 1.14$ following convention. The resulting temperature–velocity relation is given as

(2.3)\begin{equation} \tilde{T} = \tilde{T}_w + s\,Pr\, (\tilde{T}_r-\tilde{T}_w)\, \frac{\tilde{U}}{\tilde{U}_e} \left(1 - \frac{\tilde{U}}{\tilde{U}_e} \right) + \left( \frac{\tilde{U}}{\tilde{U}_e} \right)^2 ( \tilde{T}_e-\tilde{T}_w ), \end{equation}

with the recovery temperature $\tilde {T}_r = \tilde {T}_e + r \tilde {U}_e^2/(2 C_p)$. This relation has been validated across a wide range of channel flows, pipe flows and boundary layers with and without heat transfer (Zhang et al. 2014; Zhang, Duan & Choudhari 2018; Modesti & Pirozzoli 2019; Volpiani, Bernardini & Larsson 2020a; Fu et al. 2021). Specifically, this relation is derived by Zhang et al. (2014) through defining the generalized recovery temperature $\tilde {T}_{r_g} = \tilde {T} + r_g \tilde {U}^2/(2 C_p)$. Then it is assumed that $\tilde {T}_{r_g} = \tilde {T}_w + U_s \tilde {U}/C_p$, where $U_s$ is a constant velocity scale. Equivalently, the assumption can be reinterpreted that $\tilde {T}$ can be represented approximately as a second-order Taylor expansion in terms of powers of $\tilde {U}$, i.e.

(2.4)\begin{equation} \tilde{T} = b_0 + b_1 \tilde{U} + b_2 \tilde{U}^2/2, \end{equation}

where the no-slip condition implies $b_0 = \tilde {T}_w$, $b_1 = (\mathrm {d}\tilde {T} / \mathrm {d}\tilde {U})|_w$. The algebraic relation of Zhang et al. (2014) can be recovered if $b_2$ is specified by evaluating the expression at the boundary-layer edge $\tilde {T}_e=\tilde {T}|_{\tilde {U}_e}$ and $b_1$ is determined using the Reynolds analogy. However, in this work, we use the matching data (denoted with subscript $m$) $\tilde {T}_m=\tilde {T}|_{\tilde {U}_m}$ to set $b_2$, such that the exact value at the matching location can be enforced. The final temperature–velocity relation is

(2.5)\begin{equation} \tilde{T} = \tilde{T}_w + s\,Pr\, (\tilde{T}_r-\tilde{T}_w)\, \frac{\tilde{U}}{\tilde{U}_e} \left(1 - \frac{\tilde{U}}{\tilde{U}_m} \right) + \left( \frac{\tilde{U}}{\tilde{U}_m} \right)^2 ( \tilde{T}_m-\tilde{T}_w ). \end{equation}

Note that one consequence of this relation is that the wall heat flux and wall shear stress are linked algebraically by the Reynolds analogy factor, where the heat flux is defined as $q_w = s \tau _w C_p (\tilde {T}_w-\tilde {T}_r)/\tilde {U}_e$.

2.3. Implementation details

Like the classical model ((1.1)–(1.3)), the present model requires a matching temperature, velocity and density, an equation of state (the ideal gas law is used in this work and the thin boundary-layer assumption implies that the pressure is constant), and a viscosity law (either a power law or Sutherland's law, depending on the relevant reference data). In addition, the present model requires as input the velocity and temperature at the boundary-layer edge (computed using the method of Griffin, Fu & Moin 2021a) for deploying the algebraic temperature–velocity relation (2.5) due to its dependence on the recovery temperature and edge velocity. To solve the nonlinear system, the following approach is used. The incompressible ODE (1.1) with constant properties is integrated once analytically, rearranged for ${\rm d}\tilde {U}/{{\rm d} y}$, and substituted into the inverse velocity transformation (2.2) as $S$. This equation (initial-value problem with an initial guess for the wall shear stress) is solved via the shooting method, where, at each integration step, a sub-iteration determines the velocity increment that is consistent with the temperature–velocity relation (2.5) and the resulting density and viscosity at that location.

The implementation of the present model is available at the link provided in the data availability section at the end of this paper. This implementation was first developed by Griffin, Fu & Moin (2021b) to compute temperature and velocity profiles for estimating grid-point requirements in compressible flows, and this paper serves as the comprehensive documentation and the further development of the underlying inverse method for the WMLES approach for the first time. Intermediate developments were presented in Griffin, Fu & Moin (2022b), and initial results were reported in Griffin, Fu & Moin (2022a) and Griffin (2022). Kumar & Larsson (2022) used a similar procedure but with a data-driven velocity transformation (Volpiani et al. 2020b). Chen et al. (2023) and Song, Zhang & Xia (2023) approximate the mean profiles of channel flows by considering two velocity transformations (Trettel & Larsson 2016; Griffin et al. 2021c) and employing the central mean temperature scaling (Song et al. 2022).

3. A priori results

The present and classical wall models are first evaluated via a priori analysis. That is, the matching data are taken from direct numerical simulations (DNS) at wall-normal distance $y_m=0.3\delta$. The wall model estimates the velocity and temperature profiles, as well as the wall shear stress and wall heat flux. The predicted velocity and temperature profiles are shown in figures 1 and 2 for four channel flows with various Mach and Reynolds number conditions, figure 3 for two pipe flows at different Reynolds numbers, and figure 4 for two boundary layers, one with a heated and one with a cooled wall boundary condition. The bulk Mach number is defined as $M_b = U_b/\sqrt {\gamma R \tilde {T}_w}$, where $\gamma$ is the ratio of specific heats, and $R$ is the gas constant. The bulk Reynolds number is defined as $Re_b = \rho _b U_b \delta / \bar {\mu }_w$, where the bulk density is defined as $\rho _b = \iint _A \bar {\rho }\,{\rm d}A/A$, and the bulk velocity is defined as $U_b = \iint _A \tilde {U}\,{\rm d}A/A$, with $A$ the cross-sectional area of the domain. Reference DNS data are provided by Modesti & Pirozzoli (2019), Trettel & Larsson (2016), Zhang et al. (2018) and Volpiani et al. (2020a). For all cases, the profiles predicted by the present model agree with the DNS profiles significantly better than the classical model. Note that the velocities are non-dimensionalized by the predicted friction velocity, so the obtained profiles do not necessarily pass through the matching data if the predicted wall stress is inaccurate.

Figure 1. A priori wall-modelled profiles of (a,c) velocity and (b,d) temperature are plotted versus the wall-normal coordinate. Results are for supersonic channel flows with (a,b) $Re_\tau ^* = 410$, $M_b = 1.7$ and ${-B_q=0.053}$, and (c,d) $Re_\tau ^* = 590$, $M_b = 1.7$ and $-B_q=0.049$.

Figure 2. A priori wall-modelled profiles of (a,c) velocity and (b,d) temperature are plotted versus the wall-normal coordinate. Results are for supersonic channel flows with (a,b) $Re_\tau ^* = 590$, $M_b = 3.0$ and ${-B_q=0.12}$, and (c,d) $Re_\tau ^* = 200$, $M_b = 4.0$ and $-B_q=0.19$.

Figure 3. A priori wall-modelled profiles of (a,c) velocity and (b,d) temperature are plotted versus the wall-normal coordinate. Results are for supersonic pipe flows with (a,b) $Re_\tau ^* = 333.5$, $M_b = 1.500$ and $-B_q=0.047$, and (c,d) $Re_\tau ^* = 668.8$, $M_b = 1.500$ and $-B_q=0.044$.

Figure 4. A priori wall-modelled profiles of (a,c) velocity and (b,d) temperature are plotted versus the wall-normal coordinate. Results are for hypersonic diabatic (isothermal) boundary layers with (a,b) ${Re_\tau ^* = 5677}$, $M_e = 11.46$ and $-B_q=0.19$ (cooled wall), and (c,d) $Re_\tau ^* = 2328$, $M_e = 4.327$ and $-B_q=-0.039$ (heated wall).

Next, the model performance is evaluated with a wide range of DNS data from 48 different simulations. The errors in the modelled wall stress and heat flux predictions are reported for each case with $y_m=0.3\delta$. The relative error in the wall stress prediction $\epsilon _{\tau _w}$ is defined as

(3.1)\begin{equation} \epsilon_{\tau_w} = \frac{\tau_{w,{model}} - \tau_{w,{DNS}}}{\tau_{w,{DNS}}} \times 100\,\%. \end{equation}

The non-dimensional wall heat flux is defined as $B_q = q_{w}/(C_p \tilde {T}_w \bar {\rho }_w u_{\tau })$, and the relative error in the wall heat flux is defined as

(3.2)\begin{equation} \epsilon_{q_w} = \frac{ q_{w,{model}} - q_{w,{DNS}} }{q_{w,{DNS}}} \times 100\,\%. \end{equation}

Here, $\epsilon _{q_w}$ is not reported for adiabatic boundary-layer data because it is undefined, and both models predict negligible heat transfer for these data. The data considered include the compressible channel flow simulations of Modesti & Pirozzoli (2016), Trettel & Larsson (2016) and Yao & Hussain (2020), the pipe flow simulations of Modesti & Pirozzoli (2019), and the diabatic supersonic and hypersonic boundary layers of Zhang et al. (2018), Volpiani et al. (2018, 2020a) and Fu et al. (2019). The cases have edge Mach numbers in the range 0.77–11, and semi-local friction Reynolds numbers in the range 170–5700. Only the cases with $Re_\tau ^* > 150$ are analysed because lower Reynolds numbers can exhibit strong Reynolds number effects (Modesti & Pirozzoli 2016) and are not the target of this study. The error measures are shown in figure 5. The present model generates significantly less modelling error than the classical model, with the greatest error reduction when the non-dimensional heat transfer is the highest.

Figure 5. A priori modelling errors of (a) the wall shear stress $\tau _w$, and (b) the wall heat flux $q_w$, versus the heat transfer coefficient $B_q$. The model matching data are taken from DNS of various channel and pipe flows (squares) and diabatic boundary layers (circles).

To distinguish the effects of Reynolds number and compressibility, we explore the effect of using Reynolds-number-dependent coefficients for the underlying incompressible law of the wall. Specifically, rather than letting the von Kármán constant $\kappa$ and the damping coefficient $A^+$ be fixed values of 4.1 and 17, respectively, we recalibrate these values using incompressible reference data at various Reynolds numbers. We employ the DNS data from five incompressible turbulent channel flows (Lee & Moser 2015) with friction Reynolds numbers $Re_\tau = u_\tau \delta / \nu _w = \{182, 543, 1000, 1990, 5190\}$, and fit the least squares optimal values of $\kappa = \{0.400, 0.408, 0.400, 0.391, 0.391\}$ and $A^+ = \{18.2, 17.4, 17.0, 16.5, 16.5\}$. Linear interpolation and constant extrapolation of the optimal values are used to define $\kappa$ and $A^+$ for all Reynolds numbers. The inverse velocity transformation uses the semi-local wall-normal coordinate $y^*$, so the incompressible data should be interpreted as a function of $Re_\tau ^*$ rather than $Re_\tau$. A priori analysis is performed as before using compressible DNS data, but with the optimal coefficients selected according to the $Re_\tau ^*$ observed in the compressible DNS. In figures 6(a,b), for the case of a turbulent channel flow with $Re_\tau ^* = 190$ and $M_b = 1.7$, there is a modest improvement from using the Reynolds-number-dependent coefficients for the incompressible model. This suggests that at low Reynolds numbers, the deviation of DNS data for the incompressible constant-property velocity profile from the nominal law of the wall is of the same order as the deviation of the constant coefficient model and compressible DNS velocity profile. However, there is not a complete collapse of the model with Reynolds-number-dependent coefficients with the compressible DNS. This is likely attributed to the documented error in the compressible velocity transformation at $Re_\tau ^* \lesssim 200$ (Griffin et al. 2021c). In figures 6(c,d), the case of a turbulent channel flow with $Re_\tau ^* = 590$ and $M_b = 1.7$ is considered. The Reynolds number is high enough that the optimal and constant coefficients are similar; thus the performance of the present model with either set of coefficients is similar. Overall, there is no significant sensitivity to tuning the coefficients, so, for simplicity, we use the constant coefficients $\kappa =0.41$ and $A^+=17$ for the remainder of this paper.

Figure 6. A priori wall-modelled profiles of (a,c) velocity and (b,d) temperature are plotted versus the wall-normal coordinate. Results are for supersonic channel flows with (a,b) $Re_\tau ^* = 190$, $M_b = 1.7$ and ${-B_q=0.057}$, and (c,d) $Re_\tau ^* = 590$, $M_b = 1.7$ and $-B_q=0.049$.

Two more recently developed compressible wall models are considered. The first is developed by Yang & Lv (2018); they show that the damping function in the classical model (1.2) is consistent with the velocity transformation of Van Driest (1951), which has been shown to be less accurate in channel flows than the velocity transformation of Trettel & Larsson (2016). Therefore, Yang & Lv (2018) rewrite the damping function in terms of $y^*$, and show that this makes the model consistent with the Trettel–Larsson transformation. The second additional model considered is proposed by Chen et al. (2022), which also uses the semi-local damping function and further replaces the constant turbulent Prandtl number assumption of the classical model with an explicit function of $y^*$. In figure 7, these two additional wall models are compared with the classical and present wall models. Figures 7(a–d) indicate that all models are performing well in the channel flows except for the classical model. This behaviour is explained by the behaviour of the underlying velocity transformations. The models of Yang & Lv (2018) and Chen et al. (2022) use the Trettel–Larsson transformation, and the present model uses the total-stress-based transformation (Griffin et al. 2021c). Both of these transformations are well established to outperform the van Driest transformation (used by the classical model) in channel flows. In figures 7(e, f) and 7(g,h), the models are applied to boundary layers with cooled and heated walls, respectively. For both cases, the classical model is the least accurate, likely due to the inaccuracy of the van Driest transformation for boundary layers with strong heat transfer (Griffin et al. 2021c), as the velocity transformation is the only difference between the classical model and that of Yang & Lv (2018). Also for both cases, the models that use semi-local damping (Yang & Lv 2018; Chen et al. 2022) perform almost identically, suggesting limited sensitivity in these flows to the change in turbulent Prandtl number model proposed by Chen et al. (2022). For the heated boundary layer, the present model slightly improves the prediction of the temperature peak and the log slope of the velocity compared to the semi-local damping models. For the cooled boundary layer, there is a more substantial improvement from the present model for the log slope of the velocity, but the temperature profiles are only slightly improved. These improvements of the present model over the semi-local damping models are consistent with the improvements of the total-stress-based transformation over the Trettel–Larsson transformation for boundary layers with strong heat transfer.

Figure 7. A priori wall-modelled profiles of (a,c,e,g) velocity and (b,d, f,h) temperature are plotted versus the wall-normal coordinate for four wall models: (a–d) corresponding to supersonic channel flows presented in figure 1; (e–h) corresponding to hypersonic diabatic boundary layers presented in figure 4.

4. A posteriori WMLES results

In this section, several WMLES simulations are conducted using charLES, a high-fidelity compressible finite-volume code (Brès et al. 2018). The numerical method consists of a low-dissipation, approximately entropy-preserving scheme, which utilizes artificial bulk viscosity to capture the solution discontinuities. Additional details about the solver and a summary of validation campaigns are available in Fu et al. (2021, 2022).

The WMLES conducted herein are compressible turbulent channel flows driven with uniform volumetric momentum and energy source terms to achieve the same bulk Mach number $M_b$ and bulk Reynolds number $Re_b$ conditions of the DNS simulations of Trettel & Larsson (2016) as summarized in table 1.

Table 1. Non-dimensional flow parameters for the nine compressible turbulent channel flow cases considered for a posteriori testing within the WMLES framework.

The cases are run on a domain of size $({\rm \pi} \times 2 \times {\rm \pi}\sqrt {3}/4)\delta$ with periodic boundary conditions in the streamwise (first) and spanwise (third) dimensions. The mean profiles and fluxes were insensitive to doubling of the streamwise and spanwise domain sizes. Consistent with the DNS simulations, the viscosity is described by $\mu /\mu _{ref}=(T/T_{ref})^{0.75}$ and $Pr = 0.7$. All cases are initialized from a uniform solution with the target bulk Mach number and Reynolds number, and zero velocity in the wall-normal and spanwise directions. The simulations are allowed to transition from laminar to turbulent states naturally, and are run for $\sim 500$ eddy turnover times $\delta /u_\tau$. To challenge the wall model and isolate the effect of near-wall numerical errors (Kawai & Larsson 2012), the wall model matching location is placed at $y_m=0.3\delta$, and a coarse grid of 12 points per half-channel height is used for all simulations unless indicated otherwise. The computational cost of the present model is similar to that of the classical model. The present model varies between being 7 % faster and 32 % slower, depending on the Reynolds number, matching location and Mach number. No effort was made to optimize the performance of the present model, so these numbers are just meant to indicate that the approximate cost of the model is similar in the cases tested. In general, modest differences in the cost of a wall model can be amortized efficiently over parallel processors via load balancing that assigns fewer control volumes to processors that contain more boundary faces, but this is not used in the present study.

The velocity and temperature profiles from WMLES are shown in figures 8 and 9 for turbulent channel flows at four combinations of Reynolds and Mach numbers. In all cases, the present model is significantly more accurate than the classical model for the prediction of velocity and temperature with respect to the reference DNS solutions. For these cases and the others listed in table 1, the errors in the predictions of the wall shear stress and the wall heat flux are shown in figure 10. The wall model is based on the inversion of the total-stress-based velocity transformation (Griffin et al. 2021c), and that was observed to have the greatest improvement over classical approaches in cases with strong heat transfer. This explains why the errors from the classical wall model grow significantly with the strong heat transfer, but the errors from the present model are rather small and do not vary with heat flux.

Figure 8. (a,c) Velocity and (b,d) temperature profiles from WMLES with the classical (blue) and present (red) wall models. A channel flow with $M_b=1.7$, $Re_\tau ^*=410$ and $-B_q = 0.053$ is shown in (a,b), and one with $M_b=1.7$, $Re_\tau ^*=590$ and $-B_q = 0.049$ is shown in (c,d). Within the WMLES framework, the outer solutions are computed by the LES PDE solver, while the inner solutions are computed by the two wall models. These solutions coincide at the LES matching point nearest to $y_m=0.3\delta$, which is indicated with the dashed and dotted lines for the classical and present models, respectively.

Figure 9. (a,c) Velocity and (b,d) temperature profiles from WMLES with the classical (blue) and present (red) wall models. A channel flow with $M_b=3.0$, $Re_\tau ^*=590$ and $-B_q=0.12$ is shown in (a,b), and one with $M_b=4.0$, $Re_\tau ^*=200$ and $-B_q=0.19$ is shown in (c,d). Within the WMLES framework, the outer solutions are computed by the LES PDE solver, while the inner solutions are computed by the two wall models. These solutions coincide at the LES matching point nearest to $y_m=0.3\delta$, which is indicated with the dashed and dotted lines for the classical and present models, respectively.

Figure 10. The WMLES a posteriori modelling errors for (a) the wall shear stress $\tau _w$ and (b) the wall heat flux $q_w$, versus the non-dimensional heat flux $B_q$. WMLES is conducted using the classical (blue) and present (red) wall models for turbulent channel flows at the nine operating conditions listed in table 1.

The primary quantities of interest for WMLES are the predictions of the mean profiles and fluxes. The fluctuating parts of LES solutions are not expected to agree exactly with DNS results unless the WMLES are conducted with DNS-like resolution, which is impractical. Nevertheless, the effect of wall models on the fluctuating part of the LES solution is presented for comparison between the present and classical models. Figures 11 and 12 include profiles of the LES resolved turbulent Mach number $M_t=u''/\sqrt {\gamma R \tilde {T}}$ and the LES temperature fluctuations $T''$, where $''$ denotes the Favre fluctuation, with $({\cdot })'' = ({\cdot }) - \widetilde {({\cdot })}$. There is an improvement in the predictions of the fluctuating statistics by the present model compared to those by the classical model. An accurate prediction of second-order statistics is unlikely without an accurate prediction of mean statistics. Thus the improved second-order statistics of the present model are likely a consequence of its improved mean statistics compared to those of the classical model (see figures 8 and 9). However, correct prediction of the mean field is not sufficient for the accurate prediction of second-order statistics in LES. In fact, the fluctuations in the LES results are generally over-predicted compared to the DNS data. The over-prediction may be due in part to the wall-blocking effect of stress-based wall model (Bae et al. 2018). Given the coarse resolution of twelve points across the channel half-height, numerical errors and subgrid-scale model errors are certainly contributing. The subgrid-scale model has not been adapted for compressibility other than by accounting for variable properties (Moin et al. 1991). The turbulent Mach numbers are of the order of 0.3, which is sufficiently high that modelling for dilatational dissipation is a promising path to further improvements of the fluctuating statistics in the volume of the LES domain. Such research may be pursued independently of the current study focusing on wall modelling and the prediction of mean profiles and fluxes.

Figure 11. The (a,c) LES turbulent Mach number $M_t$ and (b,d) LES temperature fluctuation $T''$ profiles from WMLES with the classical (blue) and present (red) wall models. A channel flow with $M_b=1.7$, $Re_\tau ^*=410$ and $-B_q = 0.053$ is shown in (a,b), and one with $M_b=1.7$, $Re_\tau ^*=590$ and $-B_q = 0.049$ is shown in (c,d). The symbols represent the outer solutions computed by the LES PDE solver.

Figure 12. The (a,c) LES turbulent Mach number $M_t$ and (b,d) LES temperature fluctuation $T''$ profiles from WMLES with the classical (blue) and present (red) wall models. A channel flow with $M_b=3.0$, $Re_\tau ^*=590$ and $-B_q = 0.12$ is shown in (a,b), and one with $M_b=4.0$, $Re_\tau ^*=200$ and $-B_q = 0.19$ is shown in (c,d). The symbols represent the outer solutions computed by the LES PDE solver.

4.1. Sensitivity to numerical resolution and the matching location

In WMLES, the wall model exchanges data with the outer LES solver at the matching location. The modelling error in the inner wall-modelled equations may grow as the matching distance increases, which motivates placing the matching location near the wall. On the other hand, the matching location should be far enough from the wall in terms of the LES mesh resolution so that the LES solver can resolve the large scales of turbulence at the height of the matching location. Otherwise, numerical errors may contaminate the matching data that are provided as input to the wall model. Kawai & Larsson (2012) demonstrate this trade-off and how LES numerical errors contaminate the wall-modelled solution if the matching distance is of the order of the wall-normal grid resolution. The optimal matching distance will depend on the accuracy of a specific LES solver, but a typical choice is $y_m \ge 3\varDelta$ (Kawai & Larsson 2012), where $\varDelta$ is the wall-normal grid spacing near the wall.

To evaluate the convergence and sensitivity of the presently proposed wall model, two types of mesh convergence studies are considered. In the first study, the matching location is held fixed at $y_m=0.3\delta$, which corresponds in semi-local units to $y_m^*=186$ and $y_m^*=237$ for the present model and classical model cases across all resolutions. For the case with $M_b=3.0$ and $Re_\tau =1800$, the numerical resolution of the WMLES is varied. In figure 13, the WMLES solutions are shown for three LES resolutions with 9, 18 and 36 grid points across the channel half-height. The uniform hexagonally close-packed mesh topology with global refinement is employed, resulting in three meshes with $2.0\times 10^4$, $1.6\times 10^5$ and $1.3\times 10^6$ control volumes, respectively. (Note that the reference DNS uses as many as $6.4\times 10^8$ control volumes.) In this study, the LES numerical errors at the matching location are expected to diminish as the resolution is refined, but modelling errors from using the wall model over the domain $y\in [0,0.3\delta ]$ are not expected to change with resolution. For this reason, the classical model shows a large error in the log intercept of the velocity profile that is persistent with refinement and consistent with a priori analysis in figure 2(a). For the finest resolution with the present model, the grid point nearest to the wall exhibits an error that is persistent with refinement, which is consistent with the observations of (Kawai & Larsson 2012) and does not affect the accuracy of the simulation since the inner solution is applicable for $y< y_m$. For both the present and classical models, the results are only weakly dependent on the grid resolution. This suggests that the leading source of error for the simulations with the classical wall model is in fact the wall model rather than the numerical or subgrid-scale modelling errors, even on the coarsest simulation with 9 grid points per channel half-height.

Figure 13. A posteriori mesh sensitivity study of a channel flow at $M_b=3.0$ and $Re_\tau ^*=590$ with the matching location fixed at $y_m=0.3\delta$ for all cases, as indicated by the vertical dashed lines. The colours indicate the numerical resolution $\varDelta$. The outer WMLES solutions and the inner wall-modelled velocity profiles are indicated with symbols and dotted curves, respectively, for (a) the present wall model and (b) the classical wall model.

In the second grid convergence study, the models are tested in the way that WMLES is typically used in practice. That is, the matching distance is moved towards the wall as the grid is refined. In this study, two channel flows with different Reynolds number conditions are considered for three LES resolutions with 12, 24 and 48 grid points across the channel half-height. The matching locations are $y_m= 0.3\delta$, $0.15\delta$ and $0.075\delta$, respectively, which corresponds to $y_m = 4 \varDelta$ for all cases, thus the effect of near-wall LES numerical errors is expected to be minor (Kawai & Larsson 2012). In figure 14, the convergence study is performed for $M_b=3.0$ and $Re_\tau ^*=590$, and a lower Reynolds number case with $M_b=3.0$ and $Re_\tau ^*=200$ is shown in figure 15. In both cases, the accuracy of the present model is relatively high and insensitive to mesh resolution compared to that of the classical model. For the higher Reynolds number test, the matching locations in semi-local units are always in the logarithmic region of the boundary layer. Therefore, the WMLES results are not sensitive to refinement over this range of resolutions. However, for the lower Reynolds number case, the most refined meshes lead to semi-local matching locations $y_m^*$ in the buffer region. For the classical model, because the relative error of the modelled $U^+$ versus the DNS $U^+$ is maximal in the region of the buffer layer and early log layer (compare to similar a priori results in figure 6), the convergence behaviour for the classical model is complex in this regime. In other words, as the mesh is refined, although the LES numerical errors are diminishing, the wall modelling errors for the classical model may increase or decrease depending on the matching location since the relative modelling error does not reduce monotonically with wall-normal distance. On the other hand, the outer solution of the present model is relatively accurate irrespective of the matching location because the inner wall-modelled solution agrees well with the DNS solution throughout the viscous sublayer, buffer layer and log layer (which is consistent with similar a priori results in figure 6).

Figure 14. A posteriori mesh sensitivity study of a channel flow at $M_b=3.0$ and $Re_\tau ^*=590$ with matching locations dependent on the grid resolution as $y_m = 4\varDelta$. The colours indicate semi-local matching distance $y_m^*$, which is also indicated with vertical dashed lines. The outer WMLES solutions and the inner wall-modelled velocity profiles are indicated with symbols and dotted curves, respectively, for (a) the present wall model and (b) the classical wall model.

Figure 15. A posteriori mesh sensitivity study of a channel flow at $M_b=3.0$ and $Re_\tau ^*=200$ with matching locations dependent on the grid resolution as $y_m = 4\varDelta$. The colours indicate semi-local matching distance $y_m^*$, which is also indicated with vertical dashed lines. The outer WMLES solutions and the inner wall-modelled velocity profiles are indicated with symbols and dotted curves, respectively, for (a) the present wall model and (b) the classical wall model.

5. Conclusion

In this work, a wall model is proposed for turbulent wall-bounded flows with heat transfer. The model uses an established ODE description of incompressible flow, transforms that equation to account for compressibility effects, and is closed with an algebraic temperature–velocity relation. The resulting model can estimate accurately the near-wall profiles of temperature and velocity when the matching location is in the inner layer. This model is suitable for deployment as a boundary condition for an outer LES or RANS solver, an inflow generation scheme, or the base flow for perturbation methods, possibly with the incompressible model augmented with a wake profile for the outer layer of the boundary layer. The proposed method can only be as accurate as the models on which it is based, namely, the forward velocity transformation and the algebraic temperature–velocity relation. While these models have been validated widely in channel and pipe flows, and boundary layers with moderate pressure gradients, further studies in complex flows are warranted, e.g. the developing boundary layers on a blunt body behind a curved shock.

The model is first tested a priori to verify that it can recover the boundary-layer velocity and temperature data when provided with matching data from DNS. Numerical results reveal that the model recovers accurately the targeted profiles, and the predicted wall stress and heat flux are within a few per cent of their expected values for a wide database of DNS data for high-Mach-number turbulent channel flows, pipe flows and boundary layers (48 cases, with edge Mach numbers in the range 0.77–11, and semi-local friction Reynolds numbers in the range 170–5700). The model is also tested a posteriori as a boundary condition for WMLES in turbulent channel flows with bulk Mach numbers $M_b=0.7$–$4.0$ and $Re_\tau =320$–$1800$. Especially in flows with strong heat transfer, the proposed model is substantially more accurate than the classical ODE-based near-wall model. The superior performance of the present model is due to two key differences with respect to the classical model: (1) the constant turbulent Prandtl number assumption is replaced with a more accurate algebraic temperature–velocity relation; and (2) the van Driest velocity transformation is replaced with the total-shear-stress velocity transformation (Griffin et al. 2021c).