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

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 simulation (LES), high spatial and temporal resolutions are required to accurately simulate 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-modeled (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 modeling (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 widely used as boundary conditions for under-resolved RANS simulations (e.g., Abrahamson & Brower (1988); Lien et al. (1998); Goncalves & Houdeville (2001); Parente et al. (2011)). In this context, such a model is typically referred to as a wall function. Cabot (1995); 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,

where $\overline{\mu}$, $\overline{\mu}_{t}$, and $\widetilde{U}$ are the molecular dynamic viscosity, eddy viscosity, and velocity profiles, respectively, and $y$ is the wall-normal coordinate. $\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 modeled as

where $\overline{\rho}(y)$ is the density profile. The subscript $(\cdot)_{w}$ denotes quantities evaluated at the wall. $\tau_{w}=\overline{\mu}_{w}(d\widetilde{U}/dy)_{w}$ is the wall shear stress. The superscript $(\cdot)^{+}$ denotes non-dimensionalization by the friction velocity $u_{\tau}=\sqrt{\tau_{w}/\overline{\rho}_{w}}$, $\overline{\rho}_{w}$, and the kinematic wall viscosity $\overline{\nu}_{w}=\overline{\mu}_{w}/\overline{\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 ODE in Eq. (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 (Eq. (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 PDE (outer) simulation.

For compressible flow, the RANS equation for temperature can similarly be simplified 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.,

where $\widetilde{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 Eq. (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 non-linearity 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 et al., 2022), however, Iyer & Malik (2019) reports 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 et al., 2021c).

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

There are two principal differences between the present model and the classical ODE-based wall model (Eqs. (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.

The present and classical wall models are first evaluated via a priori analysis. That is, the matching data are taken from DNS at a wall-normal distance of $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 Figure 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\widetilde{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/\overline{\mu}_{w}$, where the bulk density is defined as $\rho_{b}=\iint_{A}\overline{\rho}dA/A$ and the bulk velocity is defined as $U_{b}=\iint_{A}\widetilde{U}dA/A$, where $A$ is the cross-sectional area of the domain. Reference DNS data are provided by Modesti & Pirozzoli (2019); Trettel & Larsson (2016); Zhang et al. (2018); 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.

Next, the model performance is evaluated with a wide range of DNS data from 48 different simulations. The errors in the modeled 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

The non-dimensional wall heat flux is defined as $B_{q}=q_{w}/(C_{p}\widetilde{T}_{w}\overline{\rho}_{w}u_{\tau})$, and the relative error in the wall heat flux is defined as

$\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); Yao & Hussain (2020), the pipe flow simulations of Modesti & Pirozzoli (2019), the adiabatic supersonic and hypersonic boundary layers of Pirozzoli & Bernardini (2011); Zhang et al. (2018); Volpiani et al. (2018, 2020a), and the diabatic supersonic and hypersonic boundary layers of Zhang et al. (2018); Volpiani et al. (2018, 2020a); Fu et al. (2019). The cases have edge Mach numbers in the range of 0.77–11 and semi-local friction Reynolds numbers in the range of 170–5700. Only the cases with $Re_{\tau}^{*}>150$ are analyzed 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 modeling error than the classical model, with the greatest error reduction when the non-dimensional heat transfer is the highest.

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 Figure 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 on 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}^{*}<\sim 200$ Griffin et al. (2021c). In Figure 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 of $\kappa=0.41$ and $A^{+}=17$ for the remainder of this manuscript.

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 (Eq. (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. Figure 7(a-d) indicate that all models are performing well in the channel flows except for the classical model. This behavior is explained by the behavior 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.

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 WMLESs 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.

The cases are run on a domain of size $(\pi\times 2\times\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 otherwise indicated. 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 efficiently amortized 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 Figure 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 12. 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.

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 exactly agree with DNS results unless the WMLES is 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 10 and 11 include profiles of the LES resolved turbulent Mach number $M_{t}=u^{\prime\prime}/\sqrt{(}\gamma R\tilde{T})$ and the LES temperature fluctuations $T^{\prime\prime}$, where $(\cdot)^{\prime\prime}$ denotes the Favre fluctuation $(\cdot)^{\prime\prime}=(\cdot)-\tilde{(\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 Figure 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 on the order of 0.3, which is sufficiently high that modeling 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 modeling and the prediction of mean profiles and fluxes.

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 accurately estimate 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 widely validated 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 accurately recovers the targeted profiles well, and the predicted wall stress and heat flux are within a few percent 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 of 0.77–11 and semi-local friction Reynolds numbers in the range of 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).

Kevin Griffin acknowledges support from the National Defense Science and Engineering Graduate Fellowship, the Stanford Graduate Fellowship, the Stanford Lieberman Fellowship, and the Exascale Computing Project (Grant17-SC-20SC), a collaborative effort of two US Department of Energy organizations (Office of Science and the National Nuclear Security Administration) responsible for the planning and preparation of a capable exascale ecosystem, including software, applications, hardware, advanced system engineering, and early testbed platforms, in support of the nation’s exascale computing imperative. Lin Fu acknowledges funding from the Research Grants Council (RGC) of the Government of Hong Kong Special Administrative Region (HKSAR) with RGC/ECS Project (No. 26200222) and from the Guangdong Basic and Applied Basic Research Foundation (No. 2022A1515011779). Parviz Moin acknowledges support from NASA grant (No. NNX15AU93A). We wish to gratefully acknowledge helpful comments from Sanjeeb T. Bose.

This work was authored in part by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. The views expressed in the article do not necessarily represent the views of the DOE or the U.S. Government. The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. Government purposes.

[Declaration of interests]The authors declare that they do not have any financial or non-financial conflict of interests.

[Data availability statement]The data that support the findings of this study are available from the corresponding authors upon reasonable request. Matlab code implementing the proposed model will be available in the following public repository after the manuscript is accepted for publication: https://github.com/kevingriffin1/comp_wm

[Author ORCID] Kevin Griffin (0000-0002-0866-6224); Lin Fu (0000-0001-8979-8415)