2.1. Numerical solver and boundary conditions

The simulations presented in this study are conducted using the finite-volume compressible solver charLES, which computes the solution on arbitrary polyhedral meshes. Specifically, charLES utilizes a low-dissipation spatial discretization based on principles of discrete entropy preservation (Tadmor Reference Tadmor2003; Chandrashekar Reference Chandrashekar2013), in which the fluxes are constructed to globally conserve entropy in inviscid shock-free flows, and to conserve the kinetic energy in inviscid low-Mach-number flows. Artificial diffusivity is employed in order to suppress oscillations in the vicinity of shock waves. Conserved quantities (i.e. mass, momentum and total energy) are explicitly integrated in time using a three-stage strong-stability-preserving Runge–Kutta scheme (Gottlieb, Shu & Tadmor Reference Gottlieb, Shu and Tadmor2001). The spatial and temporal schemes converge to second order and third order with respect to the nominal mesh spacing and time step, respectively. Additional discussions regarding the solver discretization and its capabilities can be found in Lozano-Durán, Bose & Moin (Reference Lozano-Durán, Bose and Moin2020), Lakebrink et al. (Reference Lakebrink, Mani, Rolfe, Spyropoulos, Philips, Bose and Mace2019), Bres et al. (Reference Bres, Bose, Emory, Ham, Schmidt, Rigas and Colonius2018), and Lehmkuhl et al. (Reference Lehmkuhl, Park, Bose and Moin2018). A set of validation and verification exercises for charLES is provided in appendix A that includes hypersonic laminar boundary layers, evolution of small amplitude disturbances in a high-Mach-number channel flow, along with a hypersonic flow around the boundary-layer transition (BOLT) subscale vehicle geometry.

The formulation of the problem is described in Yang et al. (Reference Yang, Urzay, Bose and Moin2017b). Briefly, the charLES code integrates the conservation equations of mass, momentum and total energy. Favre-filtered versions of these equations are employed for large-eddy simulation (LES) cases, with the subgrid-scale (SGS) tensor and SGS energy flux being modelled using the constant-coefficient Vreman model (Vreman Reference Vreman2004), with model constant $0.07$, along with a constant SGS turbulent Prandtl number $Pr_{sgs}=0.90$. The conservation equations are supplemented with Sutherland's law for the dynamic viscosity under a constant molecular Prandtl number $Pr=0.72$ (with Sutherland's model constants satisfying $T_{{ref}} = T_{1}$ and $S/T_{1} = 1.69$, with $T_1$ being the temperature of the inflow free stream), the ideal gas equation of state and the assumption of calorically perfect gas with $\gamma =1.4$.

The geometry and operating conditions are explained in Schülein (Reference Schülein2014), Sandham et al. (Reference Sandham, Schülein, Wagner, Willems and Steelant2014) and Willems et al. (Reference Willems, Gülhan and Steelant2015). Specifically, air at Mach $Ma_{1}=U_1/a_1=6.0$, based on the inflow free-stream velocity $U_1$ and speed of sound $a_1$, flows over an isothermal flat plate held at temperature $T_w=4.5T_{1}$, as schematically shown in figure 2. In these conditions, in which $T_w$ is smaller than the free-stream stagnation temperature $T_{0}$ (i.e. $T_w/T_{0}=0.55$), the plate behaves as a cold one that receives heat from the flow. The resulting temperature profile in the wall-normal direction is non-monotonic, which is challenging to resolve with WMLES-like coarse grid resolution near the wall, as sketched in figure 1. A wedge held above the plate is responsible for generating the shock wave that impinges on the boundary layer. In this work four wedge angles $\alpha =5^\circ$, $6^\circ$, $7^\circ$ and $8^\circ$ are studied, while keeping all other parameters constant. However, the wedge is not explicitly included in the computational domain, and, therefore, the expansion fan generated by its trailing edge is not considered. Instead, the shock wave emanating from the leading edge of the wedge is imposed by appropriate jump boundary conditions, as described below.

The Cartesian coordinate system $\{x,y,z\}$ used for the analysis is shown in figure 2, with $x=0$ corresponding to the leading edge of the plate. At the inlet of the computational domain, the Reynolds number is $Re_{1,\delta _1^\star } =U_{1}\delta _1^\star /\nu _{1}= 6830$ based on the inflow values of the displacement thickness $\delta _1^\star$ and of the free-stream velocity $U_1$ and kinematic viscosity $\nu _1$. The Reynolds number based on the distance $x_1=46\delta _1^\star$ from the leading edge of the plate to the inlet plane is $Re_{1,x_1} =U_1 x_1/\nu _1= 314\,252$. Correspondingly, the similarity solution for compressible laminar boundary layers is imposed at the inlet. In addition, periodic boundary conditions are used in the spanwise direction, while a characteristic non-reflecting boundary condition, with reference pressure chosen equal to the free-stream pressure, is applied at the outlet at a downstream distance $x_e$ such that $(x_e-x_1)/\delta _1^\star =600$, where the Reynolds number based on the inflow free-stream conditions is $Re_{1,x_e} =U_1 x_e/\nu _1= 4\,410\,887$. Note that the dimensionless streamwise distance from the edge of the plate, $(x-x_1)/\delta _1^\star$, and the Reynolds number based on the streamwise coordinate, $Re_{1,x}=U_1x/\nu _1$, can be used interchangeably for quantifying the streamwise distance in the plots below by using the relation

(2.1)\begin{equation} Re_{1,x}=\left(\frac{x-x_1}{\delta_1^\star}\right)Re_{1,\delta_1^\star}+Re_{1,x_1}.\end{equation}

Different free-stream conditions emerge downstream of the recompression shock, denoted below by the subscript ‘$2$’, as in $U_2$, $\rho _2$, $T_2$, $a_2$ and $\nu _2$. These quantities are useful, for instance, when examining the turbulent boundary layer ensuing downstream of the interaction, and they are utilized later in the text for defining the post-interaction values of the Reynolds number $Re_{2,x}=U_2x/\nu _2$ and Mach number $Ma_2=U_2/a_2$.

For a given wedge angle $\alpha$, the shock is made to emanate downwards from the top boundary of the domain at a streamwise position $x_{s}$ such that the point of inviscid intersection between the shock and the plate is located at a streamwise distance, $x_{{imp}}$, is given by $(x_{{imp}}-x_1)/\delta _1^\star =350$ in all cases, where the Reynolds number is $Re_{1,x_{{imp}}}=U_{1}x_{{imp}}/\nu _{1}=2\,704\,752$, as indicated in figure 2. For $x>x_{s}$, an oblique flow entering the domain is prescribed at the top boundary using the Rankine–Hugoniot jump conditions for pressure, density and velocities at the corresponding shock strength determined by the wedge angle $\alpha$, while the discretized fluxes at the boundary cell faces are obtained by solving a Riemann problem with a Harten–Lax–van Leer-contact solver. The similarity solution for the compressible laminar boundary layer is imposed at the top boundary for $x<x_{s}$, including the vertical displacement velocity.

The simulations were initialized using the similarity solution for the laminar compressible boundary layer in the absence of an incident shock, and were evolved for 50 flow-through times. Cumulative statistics were calculated based on an on-the-fly analysis of the solution at every time step during six and eight flow-through times in DNS and WMLES, respectively. In the notation below, $\bar {f}$ and $\tilde {f}$ denote, respectively, Reynolds and Favre averages of $f$, whereas $f'=f-\bar {f}$ and $f^{''}=f-\tilde {f}$ are the corresponding fluctuations.

2.2. Computational grids

The dimensions of the computational domain are $600\delta _1^\star \times 75\delta _1^\star \times 45\delta _1^\star$ in the streamwise, wall-normal, and spanwise directions, respectively. The Cartesian grid used for DNS is $6000 \times 600 \times 400$ (1440 million cells) and is stretched in the wall-normal direction using a hyperbolic tangent clustering with a ratio of ${\rm \Delta} y_{top}/{\rm \Delta} y_w=10$. The resolution of the DNS grid utilized here is comparable to the grid resolution employed in other studies on spatially evolving compressible turbulent boundary layers, including Sandham et al. (Reference Sandham, Schülein, Wagner, Willems and Steelant2014), Adams (Reference Adams2000), Volpiani, Bernardini & Larsson (Reference Volpiani, Bernardini and Larsson2018), Pirozzoli, Bernardini & Grasso (Reference Pirozzoli, Bernardini and Grasso2010) and Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011). In addition, the DNS grid resolution employed here leads to reasonable agreement of statistical quantities such as the skin friction and the velocity-temperature relation with well-established correlations.

Two different uniform Cartesian meshes are used for the WMLES to study the effects of grid resolution in the main text. The baseline WMLES grid is $1024 \times 270 \times 144$ (40 million cells), whereas the coarse WMLES grid employs a coarser resolution in the wall-normal direction and is $1024 \times 192 \times 144$ (28 million cells) in order to assess the effects of varying the matching location between the wall model and the outer LES. The near-wall resolution in viscous units is listed in table 1 for all cases. In the baseline WMLES cases, the boundary layer was resolved with five points across the inlet plane (four points in the coarse WMLES), 27 points across the outlet plane (nine points in the coarse WMLES), and 11 points across the wall-normal plane intersecting the streamwise location of maximum wall heat flux (seven points in the coarse WMLES) or, equivalently, at the streamwise location of maximum Stanton number $St$, the latter being formally defined in § 3.

Table 1. Minimum grid spacing near the wall in viscous units $\overline {\nu _w}/u_\tau$ at the outlet of the computational domain. In this notation, $\overline {\nu _w}$ is the time- and spanwise-averaged kinematic viscosity at the wall and $u_\tau =\sqrt {\overline {\tau _w}/\overline {\rho _w}}$ is the friction velocity based on time- and spanwise-averaged values of the wall shear stress $\overline {\tau _w}$ and density at the wall $\overline {\rho _w}$.

2.3. Wall-model parameters

In the WMLES cases the equilibrium wall model described in appendix B (see also Kawai & Larsson Reference Kawai and Larsson2012; Yang et al. Reference Yang, Urzay, Bose and Moin2017b) is utilized within a wall-modelled layer adjacent to the wall. Briefly, the equilibrium wall model consists of localized, RANS-like, steady one-dimensional versions of the wall-parallel momentum equation and the stagnation energy equation for a calorically perfect gas, with eddy-viscosity closures for the turbulent transport of momentum and energy, the latter relying on the assumption of a constant turbulent Prandtl number of $0.90$. A van Driest damping function with constant $A^{+}=17$ is employed to exponentially suppress the eddy viscosity for $y^{+}\lesssim A^{+}$ in favour of the molecular viscosity. Friction scaling is employed for the van Driest damping function, since mean density variations introduced by semi-local scaling have little effect because of the moderate wall-cooling levels utilized here. In addition, the ideal gas equation of state is utilized in the wall model to relate the density $\rho$ with the temperature $T$, in such a way that the pressure across the wall-modelled layer remains equal to the pressure at the matching location $y=h_{wm}$.

The equations of the wall model are subject to non-slip and isothermal ($T=T_w$) boundary conditions at the wall, and to the instantaneous filtered values of the wall-parallel velocity, temperature and pressure at the matching location. The outputs of the wall model are the local values of the wall shear stress $\tau _w$ and wall heat flux $q_w$, which are employed as boundary conditions for the LES conservation equations of the bulk flow.

The thickness of the wall-modelled layer $h_{wm}$ employed in these simulations is equivalent to a single cell of the WMLES grid. Whereas Kawai & Larsson (Reference Kawai and Larsson2012) have shown that this choice may lead to a log-layer mismatch, the results in Yang, Park & Moin (Reference Yang, Park and Moin2017a) indicate that temporal filtering alleviates this problem. In this work the approach proposed by Yang et al. (Reference Yang, Park and Moin2017a) is used because of its simplicity of implementation in unstructured grid environments.

Since the wall model does not incorporate streamwise variations of any quantity, the upstream propagation of elliptic effects within the wall-modelled region – for instance, due to the shock-induced adverse pressure gradient – can only occur through the boundary conditions applied at the matching location. As shown in figure 3, for both WMLES resolutions, the Mach number $Ma_{wm}$ based on the time- and spanwise-averaged values of the streamwise velocity and local speed of sound at the matching location is everywhere less than 0.5 in the laminar portion for the case $\alpha =7^\circ$. This is also the case for the other values of the wedge angle treated here. These considerations indicate that the wall-modelled layer is fully subsonic on average, and that the resolved field near the matching location is the one supporting the propagation of elliptic effects. Note that, had the wall-modelled layer been thick enough to bear the sonic line inside, no propagation of elliptic effects close to the wall would have been accounted for in the WMLES.

Figure 3. Distribution of the WMLES local Mach number $Ma_{wm} = [\bar {\rho } (\bar {u}^2+\bar {v}^2)/(\gamma \bar {P})]^{1/2}$ at the matching location $y=h_{wm}$ for the case $\alpha =7^\circ$, based on the time- and spanwise-averaged values of the local velocity, density and pressure. The vertical dashed line denotes the inviscid shock-impingement location on the wall.

In figure 4 the matching location expressed in viscous units, $h^{+}_{wm}$, plunges at $Re_{1,x}\simeq 10^6$ for the case $\alpha =7^\circ$ because the flow separates there, and increases rapidly near the inviscid shock-impingement location $Re_{1,x_{{imp}}}$ due to the sharp rise of the skin-friction coefficient, as shown below in § 3. Whereas the time- and spanwise-averaged value of $h^{+}_{wm}$ remains everywhere around or below the damping constant $A^+$ in the baseline WMLES shown in figure 4(a), its maximum value overtakes $A^+$ by a factor of four. As a consequence, based on the averaged $h^{+}_{wm}$, it may be tempting to disregard the effects of the eddy viscosity built in the wall model in the baseline WMLES. Nonetheless, it is shown in § 3 that the baseline WMLES without eddy viscosity in the wall-model equations (i.e. $\mu _{t,wm}=0$) does not lead to satisfactory results neither in the transitional nor in the turbulent portions of the boundary layer. Despite the fact that the eddy-viscosity hypothesis is questionable in transitional scenarios, these considerations highlight its dynamical relevance in regions where local overshoots in $h^{+}_{wm}$ occur.

Figure 4. Distribution of the WMLES matching location $h^+_{wm}$ scaled in wall units for the case $\alpha =7^\circ$, minimum (dashed line), maximum (dot–dashed line), along with the time- and spanwise-averaged value (solid line). The vertical dashed line denotes the inviscid shock-impingement location on the wall. Included are the data for (a) baseline and (b) coarse WMLES cases.

In both baseline and coarse WMLES cases, the equilibrium wall model is applied everywhere along the surface of the plate, including the laminar portion of the boundary layer. Two important aspects are worth remarking with regards to this choice that are discussed in the remainder of this section.

It is shown in appendix C that the WMLES adequately captures the velocity and temperature profiles in the laminar boundary layer, which remains mostly steady and two dimensional until it becomes highly disturbed in the interaction region. The wall model performs correctly there because its conservation equations are equivalent to the steady laminar boundary-layer equations very close to the wall, where advection is negligible. This can be understood by examining the distribution of $h^{+}_{wm}$ in figure 4. The values of $h^{+}_{wm}$ in both WMLES remain much smaller than $A^+$ in the laminar region, thereby yielding negligible values of the eddy viscosity in the wall model. Since order-unity values of $h^{+}_{wm}$ in the laminar region are equivalent to very small values of $h_{wm}$ relative to the boundary-layer thickness, namely $h_{wm}/\delta ^\star _1 =O(Re_{1,x_1}^{-7/4})\ll 1$, the constant molecular stress predicted by the wall model in the first approximation for $y^{+}_{wm}/A^+\ll 1$ (i.e. see (B 1) in appendix B) is equivalent to the $y/\delta ^\star \rightarrow 0$ limit of the steady laminar boundary-layer equations in the absence of streamwise pressure gradient.

That the wall model performs correctly in the laminar portion of this flow can also be understood by noticing that the scenarios sought for transition in this study are not the classical ones in which unstable eigenmodes grow relatively slowly along the entire portion of the laminar boundary layer, eventually producing transition far downstream in a way that is rather well understood, at least for calorically perfect gases flowing over smooth flat surfaces in the absence of incident shocks (Mack Reference Mack1984). Instead, the physical processes leading to transition in the present study are spatially localized downstream of the shock on the leeward side of the separation bubble, and are triggered by the absolute instability of the separation bubble without participation of any intentional disturbances at the inflow (Hildebrand et al. Reference Hildebrand, Dwivedi, Nichols, Jovanović and Candler2018). As a result, in the present study there are lesser consequences derived from the fact that neither the coarse grid resolution in WMLES nor the equilibrium wall model itself can appropriately support the growth of eigenmodes along the lengthy laminar portion of the boundary layer upstream of the shock. The task of the wall model there is limited to providing the velocity and temperature profiles within the fully viscous wall-modelled layer.

In this work the computational cost of using WMLES was $\sim$150 times less than DNS. Specifically, typical DNS cases took 25 million core hours at Argonne's Mira supercomputer, whereas only 150 000 core hours were required on average for each WMLES case on the same machine. Furthermore, it is also shown in § 3 that non-wall-modelled LESs at the resolutions listed in table 1 provide completely wrong predictions in the transitional and fully turbulent zones of the boundary layer, which underscores the positive role of the wall model in warranting acceptable predictions.

3.1. Effects of the shock incidence angle on peak thermomechanical loads

The DNS distributions of $C_{f,1}$ and $St_1$ as a function of the streamwise Reynolds number $Re_{1,x}$ are provided in figure 5 for the wedge angles considered here. Initially all the curves collapse on the laminar correlation obtained from the similarity solution, as expected by the scaling with the pre-interaction free-stream values used in (3.1) and (3.2). The characteristic shapes of $C_{f,1}$ and $St_1$ include an early drop in the laminar zone due to boundary-layer separation and a sudden overshoot downstream of the shock-impingement region because of transition. The separation of the laminar boundary layer causes a change of sign in $C_{f,1}$ due to the flow reversal and a decrease in $St_1$ due to the resulting weaker temperature gradient at the wall. In contrast, transition to turbulence leads to large spikes in $C_{f,1}$ and $St_1$, whose magnitude increase with the wedge angle. An additional discussion of this important phenomenon is provided in § 3.2 upon examining flow structures participating in the augmentation of the local thermomechanical loads.

Figure 5. Direct numerical simulation results of (a) skin-friction coefficient and (b) Stanton number as a function of the local Reynolds number $Re_{1,x}$ and the wedge angle $\alpha$. In this figure the van Driest turbulent correlation for the skin-friction coefficient $C_{f,2}$ is calculated based on post-interaction free-stream conditions, with a virtual origin equated to the leading edge of the plate, and is then re-scaled by a factor of $\rho _2 U_2^2/\rho _1U_1^2$ obtained by the DNS solution to refer the skin-friction coefficient to the pre-interaction free stream, $C_{f,1}$. Similarly, the van Driest turbulent correlation for the Stanton number $St_2$ is calculated from $C_{f,2}$ using the Reynolds analogy factor $2St_2/C_{f,2}=Pr^{-2/3}$, and is then re-scaled by a factor of $\rho _2 U_2 (T_{aw,2}-T_w)/[\rho _1 U_1 (T_{aw,1}-T_w)]$ obtained by the DNS solution to refer the Stanton number to the pre-interaction free stream, $St_1$.

The case $\alpha =5^\circ$ behaves distinctly from the others. While overshoots are observed for higher wedge angles, the $\alpha =5^\circ$ case is characterized by a modest rise in $C_{f,1}$ and $St_1$, both of which stay far below the other cases. This is attributed to the fact that transition did not occur within the computational domain in the DNS of the $5^\circ$ case. Instead, the slight increments in $C_{f,1}$ and $St_{1}$ downstream of the shock are mainly produced by the variation of the free-stream aerothermodynamic state across the shock and its impact on $\overline {\tau _w}$ and $\overline {q_w}$, whereas the normalization used for $C_{f,1}$ and $St_1$ involves only the pre-interaction free-stream aerothermodynamic state, as indicated above.

Based on the above considerations, figure 5 suggests that the critical wedge angle for the onset of shock-induced transition is somewhere between $5^\circ$ and $6^\circ$. For the transitioning cases $\alpha =6^\circ$, $7^\circ$ and $8^\circ$, the values of $C_{f,1}$ and $St_1$ downstream of transition do not agree well with the turbulent correlation of van Driest as expected, since both $C_{f,1}$ and $St_1$ are based on the pre-interaction values of the free stream, as mentioned above. In addition, the van Driest turbulent correlation for the Stanton number makes use of the Reynolds analogy factor $Pr^{-2/3}=1.24$ traditionally used to approximate a Reynolds analogy for boundary layers with non-unity Prandtl numbers. It is shown in § 3.4 that agreement with the van Driest turbulent correlations for the skin-friction coefficient and Stanton number is obtained using the definitions (3.3) and (3.4) along with the modified Reynolds analogy factor of 1.16 proposed by Chi & Spalding (Reference Chi and Spalding1966).

Qualitative comparisons between the DNS Stanton numbers for $\alpha =6^\circ$, $7^\circ$ and $8^\circ$ in figure 5 with the experimental measurements by Currao et al. (Reference Currao, Choudhury, Gai, Neely and Buttsworth2020) for $\alpha =10^\circ$ show that (i) the minimum value of $St_1$ is in the separated region in both DNS and experiments, and (ii) a monotonic increase of $St_1$ occurs near the reattachment in both DNS and experiments, after which transition of the boundary layer takes place simultaneously with an overshoot in $St_1$. Downstream of the transition zone, $St_1$ decays in both DNS and experiments, although the decay in the latter is much more substantial because of the expansion fan emanating from the trailing edge of the wedge.

Small changes in the wedge angle have profound consequences on the flow field. In particular, the DNS results for $C_{f,1}$ and $St_1$ in figure 5 indicate that increasing the wedge angle leads to earlier boundary-layer separation, longer separation bubbles and higher overshoots of $C_{f,1}$ and $St_1$ near the shock-impingement region as a result of earlier transition. The dependency of the peak values $C_{f,1}$ and $St_1$ on the wedge angle $\alpha$ is shown in figure 6. The trend in the DNS results is nearly linear, such that a $1^\circ$ increase in $\alpha$ causes approximately a 30 % increase in the average peak thermomechanical load acting on the plate. Comparisons between DNS and WMLES predictions of peak values of $C_{f,1}$ and $St_1$ in figure 6 are deferred to § 3.3.

Figure 6. Direct numerical simulation (solid lines), baseline WMLES (dotted lines) and coarse WMLES (dashed lines) peak values of the (a) skin-friction coefficient and (b) Stanton number as a function of the wedge angle $\alpha$. Since the boundary layer did not transition in the DNS of the $\alpha =5^\circ$ case, its data points are purposely disconnected from the DNS data points corresponding to the transitioning cases $\alpha =6^\circ$, $7^\circ$ and $8^\circ$.

3.2. Flow field ensued by the incidence of the shock on the boundary layer

The separation of the laminar boundary layer upstream of the shock-impingement region is induced by the adverse pressure gradient created by the incident oblique shock wave, whose effect is communicated upstream along the subsonic flow close to the wall. The time- and spanwise-averaged profiles of static pressure on the wall showing the footprint of the shock wave are provided in figure 7 for the transitioning cases. The curves are composed of three plateaus (from left to right) that correspond, respectively, to the laminar zone, the separation bubble and the turbulent zone downstream of the recompression shock. Note that the third plateau may not be present in experiments subjected to expansion effects from the trailing edge of the wedge (Currao et al. Reference Currao, Choudhury, Gai, Neely and Buttsworth2020). In the present simulations, approximately five-, six- and seven-fold overall increase in the pressure is observed across the interaction region for the cases $\alpha =6^\circ$, $7^\circ$ and $8^\circ$, respectively, mostly in agreement with the inviscid theory. Similar agreements between DNS and the inviscid theory are observed for velocity and temperature ratios in table 2.

Figure 7. Direct numerical simulation (black solid lines), baseline WMLES (green dotted lines) and coarse WMLES (blue dot–dashed lines) results for time- and spanwise-averaged profiles of the wall pressure for wedge angles of (a) $\alpha =6^\circ$, (b) $\alpha =7^\circ$ and (c) $\alpha =8^\circ$. The black dashed lines indicate the dimensionless post-interaction static pressure $P_2/P_1$ calculated assuming inviscid flow.

Table 2. Comparison of the free-stream velocity, pressure and temperature ratios across the interaction zone for the transitioning cases.

A side view of the resulting separation bubble for the case $\alpha =7^\circ$ can be approximately identified as the dark triangle-shaped region at the foot of the incident shock in figure 8(a), where the gas heats up by the reversing flow deceleration and its density reaches small values. Despite the high temperatures of the gas in the separation bubble, the wall heat flux is relatively small in this region, since the velocity gradients involved in the recirculating flow are small in comparison with those present in the laminar and turbulent portions of the boundary layer. Experimental flow visualizations by Currao et al. (Reference Currao, Choudhury, Gai, Neely and Buttsworth2020) show flow features qualitatively similar to those revealed by the density field in figure 8(a).

Figure 8. Direct numerical simulation instantaneous contours for the case $\alpha =7^\circ$ including (a) sideways view of the density field, (b) zoomed view of the post-recompression-shock region, (c) plane view of the same region at $y^{+}=100$ based on outflow conditions, along with (d) time-averaged contours of the Stanton number $St_1$ along the wall. The streamwise locations indicated by the triangles are averaged in time and along the spanwise coordinate.

In addition to the separation bubble, figure 8(a) indicates that the structure of the flow ensuing from the interaction consists of a separation shock emanating from the point of flow reversal, an expansion fan radiated from the crest of the separation bubble as the supersonic overriding flow turns downwards around it, and a recompression shock created at the point of reattachment. The incident and separation shocks intersect along a horizontal line in the spanwise direction above the separation bubble, perpendicularly to the plane of figure 8(a). The result is a regular reflection that shifts the effective interaction region downstream by approximately $75\delta _1^\star$ with respect to the inviscid shock-impingement location $x_{{imp}}$. There, the effective incidence angle $\beta$ of the shock impinging on the boundary layer is closer to $\beta \approx 11^\circ$ than to the theoretical value $\beta =14.8^\circ$ corresponding to the weak solution of an oblique shock created by an $\alpha =7^\circ$ wedge. As a consequence, the effective incidence angle of the shock is always smaller than the theoretical one predicted by the inviscid solution unless the incident shock is sufficiently weak to prevent separation.

In the cases $\alpha =6^\circ$, $7^\circ$ and $8^\circ$, the boundary layer transitions to turbulence on the leeward side of the separation bubble, shortly downstream of the time- and spanwise-averaged streamwise coordinate for reattachment. A zoomed side view of this region is provided in figure 8(b). Dynamic visualizations of the flow in this region show a persistent flapping motion of the shear layer formed between the low-speed recirculating flow within the separation bubble and the high-speed flow above. This flapping motion, in conjunction with early streaks generated shortly upstream of the reattachment point, lead to the onset of broadband turbulence at the same location where the maximum value of the Stanton number occurs, as observed in figure 8(c) and further discussed below.

The DNS distribution of the time-averaged Stanton number shown in figure 8(d) for the case $\alpha =7^\circ$ suggests the presence of quasi-stationary streaky thermal footprints of the flow onto the wall near the transition region. These structures have a spanwise wavelength of approximately $5\delta _1^\star$. These structures do not vanish by increasing the averaging time interval, as corroborated by similar streaky thermal patterns observed experimentally using infrared thermography by Currao et al. (Reference Currao, Choudhury, Gai, Neely and Buttsworth2020). It should be noted that the variation of the time-averaged Stanton number along the spanwise direction is expected in the transitional region since this is the signature of the underlying instability mechanism, which is characterized by a non-zero spanwise wavenumber along with a purely exponential growth in time at each point in space (Hildebrand et al. Reference Hildebrand, Dwivedi, Nichols, Jovanović and Candler2018).

The spanwise- and time-averaged velocity and temperature profiles in the transitional region at station $(x-x_1)/\delta ^\star _1=380$ within the separation bubble are indicated by the solid lines in figure 9. The overall flow overriding the separation bubble corresponds to an inflectional shear layer. The temperature attains a maximum at $y/\delta ^\star _1=2.15$ and attenuates towards the wall due to the cold-wall boundary condition.

Figure 9. Direct numerical simulation results for the case $\alpha =7^\circ$ including time- and spanwise-averaged profiles of (a) streamwise velocity and (b) temperature. The profiles are extracted at the streamwise locations $(x-x_1)/\delta ^\star _1=380$ (solid lines; location within the separation bubble) and $(x-x_1)/\delta ^\star _1=445$ (dashed lines; location near peak heating). In (b) a peak in the temperature profile at the station $(x-x_1)/\delta ^\star _1=445$ develops very close to the wall, as shown in the inset.

The lack of monotonicity in the temperature profile in figure 9(b) in the transitional region has important consequences on the cross-correlations between velocity and temperature fluctuations in the boundary layer. To visualize this, consider the time-averaged spatial fluctuations of the streamwise and wall-normal velocities in the cross-stream plane shown in figure 10(a). Similarly to the stationary spanwise structures of the Stanton number observed in figure 8(d), the spatial inhomogeneity of the velocity fluctuations in the spanwise direction is stationary, since both are signatures of the underlying instability mechanism. Four sets of high- and low-speed streaks are observed in figure 10(a), with maximum magnitudes in the region of strong shear (${3\lesssim y/\delta ^\star _1\lesssim 5}$). The time-averaged spatial fluctuations of the streamwise and wall-normal velocities are anti-correlated along the span, where the interaction between the mean shear and streamwise vortices results in streamwise velocity streaks. The time-averaged spatial fluctuations of the temperature and streamwise velocity on the cross-stream plane are shown in figure 10(b). Two clearly distinguished regions are observed there: a first region above the wall-normal location of maximum mean temperature ($y/\delta ^\star _1=2.15$), where the fluctuations of $u$ and $T$ are anti-correlated, and a second region below the aforementioned wall-normal location, where the fluctuations of $T$ flip their sign and become positively correlated with the fluctuations of $u$. This sign change is in agreement with the lift-up effect, since positive wall-normal velocity shifts hot gas away from the wall for the portion $y/\delta ^\star _1>2.15$ of the temperature profile that has a negative gradient, whereas the opposite happens for the portion $y/\delta ^\star _1<2.15$ of the temperature profile that has a positive gradient.

Figure 10. Direct numerical simulation results for the case $\alpha =7^\circ$ including solid contours of the time-averaged fluctuations of streamwise velocity (denoted by $\langle u \rangle _t$) at $(x-x_1)/\delta ^\star _1=380$ within the separation bubble. The solid (dashed) lines indicate positive (negative) time-averaged fluctuations of the wall-normal velocity with contour spacing of 0.03 in dimensionless units (in panel (a)), as well as positive (negative) time-averaged temperature fluctuations with contour spacing of 0.05 in dimensionless units (in panel (b)).

Similar considerations as those made above also apply downstream of the separation bubble. In particular, the dashed lines in figure 9 show the time- and spanwise-averaged velocity and temperature profiles for the downstream station $(x-x_1)/\delta ^\star _1=445$ deep in the transitional zone where the Stanton number attains its maximum value. In figure 9(b) the temperature profile arrives at the wall with a positive slope because of a spike that is not visible in this vertical scale but is revealed later in § 3.4 by zoomed-up views near the wall. The corresponding time-averaged spatial fluctuations of the streamwise and wall-normal velocities on the cross-stream plane are shown in figure 11(a). The signature of the four sets of streamwise streaks that were observed upstream in figure 10(a) is still visible here. However, the structures are now distorted by harmonic interactions along the span (e.g. see $z/\delta ^\star _1\approx 5$ and 40), and the magnitude of the fluctuations is much smaller compared with the upstream values shown in figure 10(a). The time-averaged spatial fluctuations of the temperature and streamwise velocity in figure 11(b) are almost anti-correlated along the entire cross-section except for a very thin region close to the wall below the temperature peak, where the positive temperature gradient induces a change in the sign of the correlation. As a result, the non-monotonicity of the mean temperature caused by the wall coldness has a fundamental effect that leads to the breakdown of Morkovin's hypothesis below the wall-normal location of the maximum temperature. This aspect is further analysed in § 3.4.

Figure 11. Direct numerical simulation results for the case $\alpha =7^\circ$ including solid contours of the time-averaged fluctuations of streamwise velocity (denoted by $\langle u \rangle _t$) at $(x-x_1)/\delta ^\star _1=445$ near peak heating. The solid (dashed) lines indicate positive (negative) time-averaged fluctuations of the wall-normal velocity with contour spacing of 0.02 in dimensionless units (in panel (a)), as well as positive (negative) time-averaged temperature fluctuations with contour spacing of 0.05 in dimensionless units (in panel (b)).

The late stages of transition are visualized in figure 12 using instantaneous isosurfaces of the second invariant $Q$ of the velocity-gradient tensor for the case $\alpha =7^\circ$. The isosurfaces are coloured by the root-mean-square (r.m.s.) temperature and the incident shock is superimposed. The four large-scale spanwise structures discussed above are seen lingering between the reattachment and peak-heating location [$400\lesssim (x-x_1)/\delta ^\star _1\lesssim 440$]. Downstream of the peak-heating location, a breakdown into much smaller scales is observed. For the case $\alpha =7^\circ$, the boundary layer ensuing from the shock-induced transition approaches the outflow in a turbulent state approximately at a post-interaction free-stream Mach number $Ma_2=4.2$ and a post-interaction momentum-based Reynolds number $Re_{2,\theta }=2850$.

Figure 12. Direct numerical simulation instantaneous three-dimensional visualization of shock-induced transition in the case $\alpha =7^\circ$, including isosurfaces of the second invariant $Q$ of the velocity-gradient tensor coloured by the magnitude of the dimensionless root-mean-square temperature. The value of $Q$ on the isosurfaces is set to $Q=0.3 a_1^2/{{\delta ^\star _1}^2}$, with $a_1$ being the speed of sound in the free stream at the inflow. The triangles mark key locations based on time and spanwise-averaged quantities.

3.3. Comparisons between DNS and WMLES near the shock impingement

The DNS distributions of $C_{f,1}$ and $St_1$ as a function of the streamwise Reynolds number $Re_{1,x}$ are compared in figure 13 with those obtained using WMLES. In all cases, WMLES underpredicts the size of the separation bubble with respect to DNS. In addition, WMLES delays separation (i.e. the streamwise coordinate at separation predicted by WMLES is always larger than that predicted by DNS). A good agreement is however observed between DNS and WMLES in the distributions of $C_{f,1}$ and $St_1$ upstream and within the separation zone.

Figure 13. (a,c,e,g) Skin-friction coefficient and (b,d,f,h) Stanton number as a function of the local Reynolds number $Re_{1,x}$ for wedge angles (a,b) $\alpha = 5^\circ$, (c,d) $\alpha = 6^\circ$, (e,f) $\alpha = 7^\circ$ and (g,h) $\alpha = 8^\circ$. Black solid and dotted lines denote, respectively, the laminar correlation and the van Driest turbulent correlation (see figure 5 caption for details on the calculation of the turbulent correlation). Blue, yellow and red lines denote, respectively, DNS, baseline WMLES and coarse WMLES.

An additional aspect revealed by figure 13 is that WMLES prompts transition and peak heating with respect to DNS in all cases (i.e. transition starts always earlier along the streamwise coordinate in WMLES). The discrepancies between WMLES and DNS become clearly evident in the case $\alpha =5^\circ$. Specifically, the skin friction and Stanton number in WMLES rise to values significantly larger than the DNS. Closer examination of the solution in terms of the temperature contours in figure 14(a) shows that the post-interaction boundary layer undergoes transition in both baseline and coarse WMLES, whereas no transition is observed in DNS within the present computational domain. At larger wedge angles, $\alpha =6^\circ$, $7^\circ$ and $8^\circ$, both DNS and WMLES predict transition to turbulence shortly downstream of the impingement by the shock, although WMLES always does it slightly in advance with respect to DNS, as shown in figure 14(b–d).

Figure 14. Direct numerical simulation and WMLES instantaneous temperature contours on a plane parallel to the wall at $y/{\delta ^\star _1}= 1.65$ for (a) $\alpha =5^\circ$, (b) $\alpha =6^\circ$, (c) $\alpha =7^\circ$ and (d) ${\alpha =8^\circ}$. The streamwise locations indicated by triangles have been averaged in time and spanwise direction.

Apparently, the incorrect transition predicted by WMLES at $\alpha =5^\circ$, when the shock has a relatively modest effect, is caused by numerical errors, which spuriously influence the dynamics of the post-interaction boundary layer in absence of competing disturbances of physical origin (see § 1 for a discussion about the influences of inflow disturbances on transition when the shock angle is small). In contrast, at higher incidence angles, $\alpha =6^\circ$, $7^\circ$ and $8^\circ$, when the shock has a stronger effect, the WMLES prediction of transition is reasonable, albeit spatially advanced with respect to DNS. That the transition predicted by WMLES at $\alpha =6^\circ$, $7^\circ$ and $8^\circ$ at both grid resolutions is not the result of a confabulation of numerical and modelling mishaps, is evidenced, for instance, by the correct spatial trend of the transition front moving upstream as $\alpha$ increases in figure 14, or by the increase of the peak thermal load with $\alpha$ in figure 6(b) in a manner that resembles the DNS.

The distribution of the mean pressure along the wall is predicted reasonably well by WMLES, as shown in figure 7. Mismatches of 10 % are observed in the transitional region and near the separation point. In addition, both WMLES and DNS agree well with the post-interaction pressure, velocity and temperature anticipated by the inviscid solution, as observed in table 2.

The effect of the WMLES grid resolution is most significant in the transitional region. Specifically, a decrease in the WMLES grid resolution leads to shallower rises of wall pressure, $C_{f,1}$, and $St_1$ in the transitional region, along with smaller peak values of $C_{f,1}$ and $St_1$. Comparisons between DNS and WMLES peak values of $C_{f,1}$ and $St_1$ as a function of the wedge angle $\alpha$ in figure 6 indicate that improved agreement is obtained with the DNS results as the WMLES grid resolution is increased. Farther downstream of transition and peak heating, where the boundary layer becomes turbulent, only moderate differences are observed between the two resolutions, and the values of $C_{f,1}$ and $St_1$ predicted by WMLES nearly collapse on those of DNS. Overall, the baseline WMLES is closer to the DNS results in the transitional region than the coarse WMLES. Additional results are provided in appendix D, where the WMLES grid is coarsened isotropically in the three directions for $\alpha =7^\circ$.

A more detailed comparison is made between the DNS and WMLES flow fields in figure 14 by examining the instantaneous temperature contours on a plane parallel to the wall at $y/\delta ^\star _1=1.65$. In the case $\alpha =5^\circ$ shown in figure 14(a), the flow in the DNS contains organized streaks in the post-interaction region that persist downstream without undergoing breakdown. In contrast, the boundary layer in both WMLES cases involves unstable narrower streaks, and eventually transitions to turbulence. In the cases $\alpha =6^\circ$, $7^\circ$ and $8^\circ$, the discrepancies are less significant, with both DNS and WMLES leading to transition. However, in those cases, narrower streaks are observed as the grid resolution of the WMLES is increased, although these structural discrepancies do not translate into severe mismatches neither in the location of breakdown nor in the location of peak heating. The latter, indicated by red triangles in figure 14 for each case, serves as an accurate indicator of full breakdown to turbulence in both DNS and WMLES.

Further insights into the performance of the baseline WMLES are gained in figure 15 by comparing its time- and spanwise-averaged profiles of velocity and temperature near transition and peak heating with those from DNS for the case $\alpha =7^\circ$. There, the shear layer, which separates the recirculating flow and the downwards high-speed inviscid stream at the foot of the incident shock, is much closer to the wall in the WMLES, since the reattachment occurs noticeably more upstream in WMLES than in DNS (see figure 8(b) for spatial localization of the shear layer). The corresponding profiles of the r.m.s. of the fluctuations of temperature and streamwise velocity are presented in figure 16. The maximum values of those quantities occur in the shear layer in both WMLES and DNS. The turbulent heat fluxes in the streamwise and wall-normal directions, along with the Reynolds stress, are also maximized by the turbulent transport in the shear layer, as shown in figure 17. In all these profiles it is observed that WMLES provides an acceptable prediction of features such as shear-layer thickness, maximum shear, maximum values of r.m.s. fluctuations of temperature and streamwise velocity, and maximum values of turbulent heat fluxes and Reynolds stress. However, the WMLES profiles are clearly shifted in space near the shock impingement owing to different separation and reattachment locations.

Figure 15. Direct numerical simulation (solid lines) and baseline WMLES (dot–dashed lines) of time- and spanwise-averaged profiles of (a) streamwise velocity and (b) temperature at several stations along the $x$-axis for $\alpha =7^\circ$. In panel (a) the velocity is plotted as $10(\bar {u}/U_{1})+(x-x_1)/\delta _1^\star$, whereas in panel (b) the temperature is plotted as $(3/2)(\bar {T} /T_{1})+(x-x_1)/\delta _1^\star$. Also included are the solutions of the equilibrium wall model (green lines) within the wall-modelled region $y\leq h_{wm}$.

Figure 16. Direct numerical simulation (solid lines) and baseline WMLES (dot–dashed lines) of r.m.s. profiles of (a) streamwise velocity and (b) temperature at several stations along the $x$-axis for $\alpha =7^\circ$. In panel (a) the velocity is plotted as $36\sqrt {\overline {u'^2}} /U_{1}+(x-x_1)/\delta _1^\star$, whereas in panel (b) the temperature is plotted as $6\sqrt {\overline {T'^2}} /T_{1}+(x-x_1)/\delta _1^\star$.

Figure 17. Direct numerical simulation (solid lines) and baseline WMLES (dot–dashed lines) of (a) streamwise turbulent heat flux, (b) wall-normal turbulent heat flux and (c) Reynolds stress at several stations along the $x$-axis for $\alpha =7^\circ$. The data is plotted as $20\overline {\rho u''T''} /(\rho _1U_1T_1)+(x-x_1)/\delta _1^\star$ in panel (a), as $150\overline {\rho v''T''} /(\rho _1U_1T_1)+(x-x_1)/\delta _1^\star$ in panel (b), and as $1000\overline {\rho u''v''} /(\rho _1U_1^2)+(x-x_1)/\delta _1^\star$ in panel (c).

In summary, mismatches between the WMLES and DNS statistics in the transitional region are mostly dominated by the erroneous spatial advance in the reattachment predicted by WMLES. However, the root cause of this discrepancy cannot be straightforwardly isolated. This can be understood by noticing that the spatial advance in the reattachment in WMLES is coupled with the spatial delay predicted in separation, the latter engendering a separation shock at an erroneous angle. In addition, as the boundary layer approaches separation, its velocity profile becomes more contorted and critically more under-resolved by the WMLES grid. The separation shock in WMLES then intersects the main shock, which is refracted towards the boundary layer at an erroneous angle. The resulting accumulation of errors is germane to the present configuration that involves widely different, interacting flow structures communicated by an adverse pressure gradient, and is not as severely observed in WMLES predictions of less complex configurations such as shock-free flat-plate turbulent boundary layers or turbulent channel flows. These considerations highlight the closely coupled contributions of all these phenomena in setting the overall performance of WMLES in the present problem.

In order to isolate the role of the equilibrium wall model and its eddy viscosity in predicting the spatial distributions of the skin-friction coefficient and Stanton number, two tests are performed in figure 18 for the case $\alpha =7^\circ$ using the baseline WMLES grid. In the first test, the equilibrium wall model is replaced by a non-slip boundary condition. The resulting curves are denoted by the tag ‘LES without wall model’ in figure 18. Negligible changes are observed in the laminar portion upstream of the shock-impingement zone, which indicates that the wall model does not have any significant effect there. In contrast, in the post-interaction boundary layer, the skin-friction coefficient is significantly overpredicted while the Stanton number is underpredicted by a factor of two. In conclusion, without the wall model, the predictions of the transitional and turbulent portions of the boundary layer are largely degraded due to deficient physical modelling and to numerical errors enabled by the coarse LES mesh.

Figure 18. (a) Skin-friction coefficient and (b) Stanton number as a function of the local Reynolds number $Re_{1,x}$ in DNS (solid lines), baseline WMLES (thick dot–dashed lines), baseline WMLES with eddy viscosity set to zero in the wall-model equations (thin dot–dashed lines), and no-slip LES performed on the baseline WMLES grid but without the wall model (dashed lines). All results in this figure correspond to the case $\alpha =7^\circ$.

In the second test, the equilibrium wall model is activated but its eddy viscosity $\mu _{t,wm}$, defined in (B 3) in appendix B, is turned off. The resulting curves are denoted by the tag ‘WMLES without eddy viscosity in the wall-model equations’ in figure 18. By turning off the eddy viscosity, the wall model provides only the viscous continuations of the velocity and temperature profiles within the wall-modelled region. In this case, the performance of the equilibrium wall model also deteriorates significantly. In particular, despite the fact that the average matching location $h_{wm}^+$ is within the damped spatial range of the wall-normal coordinate (i.e. see figure 4), the results in figure 18 suggest that $\mu _{t,wm}$ plays an important role in the prediction of $C_f$ and $St$ not only in the turbulent boundary layer ensuing downstream of the interaction, as expected, but also in the transitional zone near the shock impingement. Note that these considerations do not imply neither that the eddy-viscosity model (B 3) is the correct one to use in the transitional zone, nor that the eddy-viscosity hypothesis is appropriate for transitional flows, but that the eddy-viscosity model (B 3) appears to have a beneficial effect on the solution compared to setting $\mu _{t,wm}=0$. This effect is particularly relevant in transitional spots of large skin friction, where the local instantaneous values of $h_{wm}^+$ can be as high as four times the van Driest damping constant, as shown in figure 4.

3.4. The turbulent boundary layer far downstream of the impingement by the shock

While §§ 3.2 and 3.3 were focused on the transitional region, this section analyses the turbulent boundary layer downstream of the interaction region for the cases $\alpha =6^\circ$, $7^\circ$ and $8^\circ$. Table 3 summarizes relevant parameters characterizing the state of the turbulent boundary layer at a representative streamwise location $(x-x_1)/\delta _1^\star =580$, on which most of the statistical analysis outlined below is based.

Table 3. Conditions corresponding to the turbulent boundary layer at $(x-x_1)/\delta _1^\star =580$ far away downstream of the recompression shock. The table includes the post-interaction Reynolds number based on the momentum thickness $Re_{2,\theta }$, the Reynolds number based on the wall-friction velocity $Re_{2,\tau }$, the post-interaction Mach number $Ma_2$ and the ratio of the wall temperature and the post-interaction free-stream temperature $T_w/T_2$.

The physical characteristics of the ensuing turbulent boundary layer are determined by the post-interaction values of the free-stream Mach number $Ma_2$, the momentum-based Reynolds number $Re_{2,\theta }$ (or the friction Reynolds number $Re_{2,\tau }$), and the ratio of the wall temperature to the free-stream temperature $T_w/T_2$. In the three cases shown in table 3, $Ma_2$ is smaller than the inflow free-stream value $Ma_1=6$ because of the net deceleration of the free stream as it crosses the incident shock and the train of shocks and expansion fans induced by the impingement. Similarly, the net heating of the free stream across the interaction zone leads to temperature ratios $T_w/T_2$ smaller than the corresponding inflow value $T_w/T_1=4.5$. In the turbulent boundary layer the overall consequences of increasing the wedge angle are a decrease in $Ma_2$, an increase in $Re_{2,\theta }$ (or $Re_{2,\tau }$), and a decrease in $T_w/T_2$.

As indicated by the third plateau of the wall pressure shown in figure 7, the turbulent boundary layer far downstream of the impingement by the shock is one under negligible mean streamwise gradient of static pressure. The wall-normal gradient of the static pressure is similarly weak, since the wall pressure in that third plateau is well described by the free-stream static pressure calculated from an inviscid interaction, as expected from the moderate values of $Ma_2$. The remainder of this section is dedicated to assessments of analogies and hypotheses traditionally developed for zero-pressure-gradient compressible boundary layers, such as Reynolds analogies, mean velocity transformations and Morkovin's hypothesis.

The three different wedge angles considered in table 3 unfold in dimensionless space as three different sets of values of $Re_{2,\theta }$, $Ma_2$ and $T_w/T_2$. Although the sensitivity of the solutions to the particular value of $Re_{2,\theta }$ is expected to be small at the relatively large values of $Re_{2,\theta }$ considered here, some of the statistics of the different metrics described below do not collapse among the three turbulent boundary layers because of additional dependencies of the solution on $Ma_2$ and $T_w/T_2$. Notable exceptions that remain relatively robust to changes of the wedge angle within the range tested here are the skin-friction coefficient $C_{f,2}$ and the Stanton number $St_2$ downstream of the recompression shock when scaled with post-interaction free-stream values, as defined in (3.3) and (3.4). This is elicited by the improved collapse of the turbulent portion of the profiles observed in figure 19, as opposed to the significant dispersion in figure 6 when pre-interaction free-stream conditions are employed instead. These considerations suggest that the effects of the variations of $\alpha$ on the wall shear stress and wall heat flux in the turbulent portion of the boundary layer can be approximately scaled out despite the different values of $Re_{2,\theta }$, $Ma_2$ and $T_w/T_2$ in each case. In figure 19(b) the agreement between the van Driest turbulent correlation for $St_2$ and the DNS solution is greatly enhanced by using the Reynolds analogy factor $2St_2/C_{f,2}=1.16$ proposed by Chi & Spalding (Reference Chi and Spalding1966) based on correlation of experimental data for turbulent boundary layers with Mach numbers less than 5 and near-adiabatic wall boundary conditions. This is in contrast to the traditional Reynolds analogy factor $2St_2/C_{f,2}=Pr^{-2/3}=1.24$ utilized in figure 5(b) for boundary layers with non-unity Prandtl numbers, which leads to significant mismatch between the van Driest turbulent correlation for $St_2$ and the DNS solution. That the Reynolds analogy factor $2St_2/C_{f,2}=1.16$ proposed by Chi & Spalding (Reference Chi and Spalding1966) is a more accurate model of the DNS results presented here can be seen in table 4. Both DNS and WMLES results settle increasingly earlier on a value of the Reynolds analogy factor as the wedge angle increases, as observed in figure 20, because the boundary layer transitions correspondingly earlier along the streamwise coordinate. The differences between the Reynolds analogy factors predicted by WMLES and DNS, and between them and the value 1.16 experimentally correlated by Chi & Spalding (Reference Chi and Spalding1966), are small and remain within a 5 % error for all the conditions tested here. However, as the wedge angle increases, the DNS results predict a slight decrease in the mean value of the Reynolds analogy factor, whereas the trend of the WMLES results is less clear. As observed in previous experimental studies collected by Cary (Reference Cary1970) and discussed in Bradshaw (Reference Bradshaw1977), increasing the wall cooling leads to a slight decrease in the Reynolds analogy factor below that proposed by Chi & Spalding (Reference Chi and Spalding1966). Although the behaviour of the DNS results observed in table 4 as the wedge angle is increased is reminiscent of an increase in wall cooling, it should be mentioned that the ratio $T_w/T_{0}$, with $T_{0}$ being the stagnation temperature, is mostly the same in both DNS and WMLES within a 0.1 % error, and is independent of the wedge angle, since the wall temperature is fixed and the stagnation temperature of the free stream is constant across the interaction region. As a result, the decrease in the Reynolds analogy factor $2St_2/C_{f,2}$ observed as the wedge angle increases in the DNS cannot be easily reconciliated with the observations made by Cary (Reference Cary1970), and may instead be attributed to the intrinsic dependency of the solution on the parameters $Re_{2,\theta }$, $Ma_{2}$ and $T_w/T_2$ listed in table 3, which differ slightly among the three cases.

Figure 19. (a) Skin-friction coefficient and (b) Stanton number as a function of the post-interaction local Reynolds number $Re_{2,x} = U_2 x/\nu _2$ for $\alpha =6^\circ$ (black lines), $\alpha =7^\circ$ (blue lines) and $\alpha =8^\circ$ (red lines), including DNS (solid lines) and the van Driest turbulent correlations (dashed lines). In this figure the van Driest turbulent correlation for the skin-friction coefficient $C_{f,2}$ is calculated based on post-interaction free-stream conditions, with a virtual origin equated to the leading edge of the plate. Similarly, the van Driest turbulent correlation for the Stanton number $St_2$ is calculated from $C_{f,2}$ using the Reynolds analogy factor $2St_2/C_{f,2}=1.16$ from Chi & Spalding (Reference Chi and Spalding1966).

Figure 20. Direct numerical simulation (blue solid lines), baseline WMLES (green dotted lines) and coarse WMLES (red dashed lines) distributions of the Reynolds analogy factor as a function of the post-interaction local Reynolds number $Re_{2,x} = U_2 x/\nu _2$ for (a) $\alpha =6^\circ$, (b) $\alpha =7^\circ$ and (c) $\alpha =8^\circ$, along with the reference value 1.16 (black dot–dashed line) proposed by Chi & Spalding (Reference Chi and Spalding1966).

Table 4. The Reynolds analogy factor, $2St_2/C_{f,2}$, for DNS, baseline WMLES and coarse WMLES, calculated at station $(x-x_1)/\delta ^\star _1= 580$ for various wedge angles.

For all three cases, figure 21 indicates that the present DNS results best match with the temperature-velocity relation proposed by Duan & Martin (Reference Duan and Martin2011), which is nonetheless based on correlation of DNS data of a different configuration involving temporally evolving turbulent boundary layers. In contrast, the Crocco–Busemann formula for $Pr=1$ (Busemann Reference Busemann1931; Crocco Reference Crocco1932), and the Walz relation that accounts for $Pr\neq 1$ (Walz Reference Walz1962, Reference Walz1966), depart from the DNS data by amounts of order 10 % and 5 %, similarly to previous observations by Zhang et al. (Reference Zhang, Bi, Hussain and She2014) and Duan et al. (Reference Duan, Beekman and Martin2010).

Figure 21. Direct numerical simulation relation between the time- and spanwise-averaged profiles of temperature and velocity at $(x-x_1)/\delta ^\star _1=580$ for a wedge angle of (a) $\alpha =6^\circ$, (b) $\alpha =7^\circ$ and (c) $\alpha =8^\circ$. Also included are the model relations from Walz (Reference Walz1962), Busemann (Reference Busemann1931), Crocco (Reference Crocco1932) and Duan & Martin (Reference Duan and Martin2011). The values used for normalization are the edge velocity $U_2$ and edge temperature $T_2$ for the corresponding angle.

In figure 21 the model for the temperature-velocity relation proposed by Duan & Martin (Reference Duan and Martin2011) requires a calibration parameter $\theta =0.8259$, which was connected later through analysis by Zhang et al. (Reference Zhang, Bi, Hussain and She2014) with the Reynolds analogy factor multiplied by the Prandtl number, namely $\theta =2St_2Pr/C_{f,2}$. Evaluation of the latter using the DNS results in table 4 indicates that $2St_2Pr/C_{f,2}$ differs from the model parameter $\theta =0.8259$ by small amounts of order $0.9\,\%$ (for $\alpha =6^\circ$), $1.7\,\%$ (for $\alpha =7^\circ$) and $1.9\,\%$ (for $\alpha =8^\circ$), thereby corroborating the analysis made by Zhang et al. (Reference Zhang, Bi, Hussain and She2014).

The comparisons between the mean temperature-velocity relations from DNS and WMLES presented in figure 22 for the case $\alpha =7^\circ$ show an encouraging agreement over the entire range of velocities. This is despite the fact that a significant portion of the momentum of the turbulent boundary layer is unresolved by the LES grid. However, the equilibrium wall model correctly captures the DNS mean temperature-velocity relations within the wall-modelled region even in the coarser WMLES case. Some understanding of the structure of the mean streamwise velocity profile can be gained by transforming it in such a way as to resemble as much as possible the mean velocity profile of an incompressible turbulent boundary layer. This is the objective of the velocity transformations shown in figure 23, which includes those proposed by van Driest (Reference van Driest1956) and Trettel & Larsson (Reference Trettel and Larsson2016), the latter being a revision of the former to account for both viscosity and density variations in boundary layers over non-adiabatic walls. Both transformations reveal the presence of viscous- and log-like layers in the transformed velocity profiles. A lack of collapse among the transformed mean velocity profiles corresponding to the three different wedge angles in figure 23 is clearly noticeable in the outer layer, where the sensitivity of the wake parameter to changes in the post-interaction Mach numbers and heating rates appears to be significant.

Figure 22. Direct numerical simulation (black solid line), baseline WMLES (red dot–dashed line) and coarse WMLES (blue dot–dashed line) relations between the time- and spanwise-averaged profiles of temperature and velocity at $(x-x_1)/\delta ^\star _1=580$ for $\alpha =7^\circ$, including the solution predicted within the wall-modelled region $y\leq h_{wm}$ (blue and red solid lines).

Figure 23. Transformed mean streamwise velocity as a function of the wall-normal coordinate for $\alpha =6^\circ$, $7^\circ$ and $8^\circ$ at $(x-x_1)/\delta ^\star _1=580$ using the transforms by (a) van Driest (Reference van Driest1956) and (b) Trettel & Larsson (Reference Trettel and Larsson2016). The figure includes DNS results for $\alpha =6^\circ$ (black solid lines), $7^\circ$ (blue solid lines) and $8^\circ$ (red solid lines), along with the incompressible profiles in the viscous sublayer (black dotted line) and log layer (dashed lines). In panel (b), $y^\star = \bar {\rho } (\overline {\tau _w}/\bar {\rho })^{1/2}y/\bar \mu$ represents a semi-locally scaled wall-normal coordinate.

Neither one of the two transformations employed in figure 23 lead to collapse of the log-layer mean velocity profile on the incompressible log law. Specifically, figure 23 indicates that, for the three angles tested here, the effective Kármán constant of the transformed mean velocity profile is smaller than the nominal Kármán constant $0.41$ of the incompressible log law.

The transformed mean velocity profiles obtained by the WMLES agree well with those of the DNS for the most part, as suggested by the comparisons provided in figure 24 for $\alpha =7^\circ$. However, discrepancies are observed in the first and second grid points of the LES grid, where the WMLES are expected to be influenced by numerical errors. The two other angles $\alpha =6^\circ$ and $8^\circ$ lead to similar conclusions and are not included here for brevity. The differences caused by coarsening the resolution of the WMLES grid in the wall-normal direction are small and do not degrade the agreement between DNS and WMLES in any significant way.

Figure 24. Transformed mean streamwise velocity as a function of the wall-normal coordinate for $\alpha =7^\circ$ at $(x-x_1)/\delta ^\star _1=580$ using the transforms by (a) van Driest (Reference van Driest1956) and (b) Trettel & Larsson (Reference Trettel and Larsson2016). The figure includes DNS (black solid lines), baseline WMLES (red dashed and solid lines), coarse WMLES (blue dashed and solid lines), along with the incompressible profiles in the viscous sublayer (black dotted line) and log layer (dashed lines). In panel (b), $y^\star = \bar {\rho } (\overline {\tau _w}/\bar {\rho })^{1/2}y/\bar \mu$ represents a semi-locally scaled wall-normal coordinate.

The good agreement shown in figure 13 between the Stanton numbers predicted by DNS and WMLES in the turbulent boundary layer downstream of the recompression shock for the transitioning cases $\alpha =6^\circ$, $7^\circ$ and $8^\circ$ must rely on the correct WMLES prediction of the mean temperature profile near the wall. This is corroborated by the comparison between the mean temperature profiles obtained from DNS and WMLES provided in figure 25(a) for $\alpha =7^\circ$. Although discrepancies of order 10 % are observed between the DNS and WMLES mean temperature profiles at wall-normal distances $y^+$ corresponding to the log and outer layers of the transformed mean velocity profile, the wall model captures correctly the mean temperature profile in the buffer zone and in the viscous sublayer. There, the temperature reaches its maximum value because of the heat generated by friction. This maximum value is not directly resolved by the LES grid but modelled successfully by the equilibrium wall model, thereby yielding a correct approximation of the magnitude and sign of the wall heat flux. The comparisons of the mean temperatures pertaining to the two other angles $\alpha =6^\circ$ and $8^\circ$ lead to similar conclusions and are not included here for brevity.

Figure 25. Direct numerical simulation (black solid lines), baseline WMLES (red dashed lines) and coarse WMLES (blue dashed lines) wall-normal profiles of (a) time- and spanwise-averaged temperature, (b) temperature/streamwise-velocity correlation coefficient, (c) temperature/wall-normal-velocity correlation coefficient and (d) turbulent Prandtl number, all profiles being obtained for $\alpha =7^\circ$ at $(x-x_1)/\delta _1^\star =580$. In panel (a) the mean temperature predicted by the equilibrium wall model within the wall-modelled region $y\leq h_{wm}$ is indicated by the blue and red solid lines. In all panels, the horizontal dashed lines indicate the wall-normal location of the maximum mean temperature.

The Morkovin hypothesis appears to provide unsatisfactory results in the present configuration. First, whereas the Morkovin hypothesis establishes perfect anticorrelation between $T''$ and $u''$ (Morkovin Reference Morkovin1962), both DNS and WMLES unisonally indicate that $T''$ and $u''$ in the present configuration are not fully anticorrelated, as shown in figure 25(b) for $\alpha =7^\circ$. Away from the wall, this non-perfect anticorrelation is explained by the approximately $10\,\%$ fluctuations observed in the stagnation temperature across the boundary layer. In addition, as anticipated in figures 10 and 11, the sign of the temperature/streamwise velocity correlation changes near the wall at the wall-normal location where the maximum of the mean temperature is attained. A similar change of sign in the temperature/wall-normal velocity correlation is also observed in figure 25(c) for $\alpha =7^\circ$ at the same location. The WMLES results provide excellent predictions for the correlations of the temperature with the streamwise and wall-normal velocities across the entirety of the resolved portion of the turbulent boundary layer.

The turbulent Prandtl number $Pr_t$ shown in figure 25(d) for $\alpha =7^\circ$ varies between 0.7 and 1.0 across the boundary layer in both DNS and WMLES, whereas a peak value of 1.5 appears to be attained at the maximum temperature location. Although not shown here for brevity, similar conclusions about the temperature/velocity correlations and the turbulent Prandtl number also hold for the other two wedge angles $\alpha =6^{\circ }$ and $8^\circ$.

That the strong Reynolds analogy (SRA), proposed by Morkovin (Reference Morkovin1962) to relate in a directly proportional way the r.m.s. values of the streamwise velocity fluctuations and the temperature fluctuations, is not a good approximation in the present configuration is shown in figure 26(a), where departures of $\sim$50 % from SRA behaviour are observed. Over the years, the SRA has been improved in different studies that account for wall heat transfer and stagnation-temperature fluctuations. For instance, Gaviglio (Reference Gaviglio1987) proposed a revised SRA (referred to as GSRA below) by assuming that the characteristic length scales of the fluctuations of temperature and velocity are similar. In a different approach, Huang et al. (Reference Huang, Coleman and Bradshaw1995) proposed another revised SRA (referred to as HSRA below) by including the local turbulent Prandtl number on the basis of a mixing-length model. Using the DNS flow fields for $\alpha =6^\circ$, $7^\circ$ and $8^\circ$, figure 26(b,c) provides an evaluation of the GSRA and HSRA expressions found in Gaviglio (Reference Gaviglio1987) and Huang et al. (Reference Huang, Coleman and Bradshaw1995), respectively, in such a way that the total validity of the corresponding relation would imply a unity value on the vertical axis across the entire boundary layer. While the classical SRA in figure 26(a) fails to reproduce the DNS data, as also observed previously by Duan et al. (Reference Duan, Beekman and Martin2010) and Zhang et al. (Reference Zhang, Bi, Hussain and She2014), the discrepancies are greatly reduced for all wedge angles by using the GSRA. However, the performance of the HSRA can be greatly enhanced by setting the Prandtl number in the HSRA relation to $Pr_t=0.9$, as shown in figure 26(d).

Figure 26. Direct numerical simulation wall-normal distributions of the (a) SRA by Morkovin (Reference Morkovin1962), (b) modified strong Reynolds analogy GSRA by Gaviglio (Reference Gaviglio1987), (c) modified strong Reynolds analogy HSRA by Huang, Coleman & Bradshaw (Reference Huang, Coleman and Bradshaw1995) based on the local turbulent Prandtl number, and (d) modified strong Reynolds analogy HSRA by Huang et al. (Reference Huang, Coleman and Bradshaw1995) based on a constant turbulent Prandtl number equal to $0.9$. The blue, red and black lines denote, respectively, the results for $\alpha =6^\circ$, $7^\circ$ and $8^\circ$ at $(x-x_1)/\delta ^\star _1=580$. The vertical dotted lines indicate unity ratios, and correspondingly, total validity of the proposed analogy. In the labels, $Ma = [\tilde {\rho } (\tilde {u}^2+\tilde {v}^2)/(\gamma \tilde {P})]^{1/2}$ denotes the local Mach number, and $\widetilde {T_{0}}$ indicates the Favre average of the local stagnation temperature $T_{0} = T[1+(\gamma -1)(\rho |\boldsymbol {v}|^2/(2\gamma P)]$ based on the modulus of the streamwise velocity vector $|\boldsymbol {v}|$. In addition, $\delta$ is the local boundary layer thickness at $(x-x_1)/\delta ^\star _1=580$ defined as the height where $\bar {u}=0.99U_2$.

The predictive capabilities of the WMLES to recreate the SRA relations is assessed in figure 27 for the representative wedge angle of $7^\circ$. Good agreement between DNS and WMLES is observed for the SRA and HSRA in figure 27(a,b). However, the WMLES deviates significantly from the DNS in the outer portion of the boundary layer when the HSRA is used with $Pr_t=0.9$, as shown in figure 27(c). These errors are commensurate with the errors incurred by the WMLES in predicting the turbulent Prandtl number calculated a posteriori from the DNS results.

Figure 27. Direct numerical simulation (black solid lines), baseline WMLES (red dot–dashed lines) and coarse WMLES (blue dashed lines) wall-normal distributions of the (a) SRA by Morkovin (Reference Morkovin1962), (b) modified strong Reynolds analogy HSRA by Huang et al. (Reference Huang, Coleman and Bradshaw1995) based on the local turbulent Prandtl number, and (c) modified strong Reynolds analogy HSRA by Huang et al. (Reference Huang, Coleman and Bradshaw1995) based on a constant turbulent Prandtl number equal to $0.9$. All profiles correspond to the case $\alpha =7^\circ$ at $(x-x_1)/\delta ^\star _1=580$. The vertical dotted lines indicate unity ratios, and correspondingly, total validity of the analogy. In the labels, $Ma = [\tilde {\rho } (\tilde {u}^2+\tilde {v}^2)/(\gamma \tilde {P})]^{1/2}$ denotes the local Mach number, and $\widetilde {T_{0}}$ indicates the Favre average of the local stagnation temperature $T_{0} = T[1+(\gamma -1)(\rho |\boldsymbol {v}|^2/(2\gamma P)]$ based on the modulus of the streamwise velocity vector $|\boldsymbol {v}|$. In addition, $\delta$ is the local boundary layer thickness at $(x-x_1)/\delta ^\star _1=580$ defined as the height where $\bar {u}=0.99U_2$.

A comparison between profiles from DNS at several streamwise stations, $(x-x_1)/\delta ^\star _1=420$, $460$ and $500$ in the transitional region, and $570$ in the fully turbulent region, is presented in figure 28 for the representative wedge angle of $7^\circ$. The correlation between the streamwise velocity fluctuation $u''$ and the temperature fluctuation $T''$ is presented in figure 28(a), the HSRA based on the local turbulent Prandtl number is shown in figure 28(b), and the turbulent Prandtl number is given in figure 28(c). While significant variations of these metrics are observed upstream deep in the transitional region, as expected for HSRA and $Pr_t$ because of their lack of clear physical meaning there, the variations among different streamwise stations tend to decrease as the turbulent portion of the boundary layer is approached, or equivalently, as the local Reynolds number $Re_{2,x}$ increases.

Figure 28. Direct numerical simulation wall-normal profiles at several streamwise stations for $\alpha =7^\circ$, including (a) temperature/streamwise-velocity correlation coefficient, (b) modified strong Reynolds analogy HSRA by Huang et al. (Reference Huang, Coleman and Bradshaw1995) based on the local turbulent Prandtl number and (c) turbulent Prandtl number. In the labels, $Ma = [\tilde {\rho } (\tilde {u}^2+\tilde {v}^2)/(\gamma \tilde {P})]^{1/2}$ denotes the local Mach number, and $\widetilde {T_{0}}$ indicates the Favre average of the local stagnation temperature $T_{0} = T[1+(\gamma -1)(\rho |\boldsymbol {v}|^2/(2\gamma P)]$ based on the modulus of the streamwise velocity vector $|\boldsymbol {v}|$. In addition, $\delta$ is the boundary layer thickness defined at $(x-x_1)/\delta ^\star _1=580$ as the height where $\bar {u}=0.99U_2$.

The streamwise and wall-normal r.m.s. velocities, normalized with semi-local inner scalings, are presented in figure 29 for $\alpha =6^{\circ }$, $7^{\circ }$ and $8^\circ$. Similarly to the transformed mean velocities in figure 23, a collapse of the curves corresponding to the three wedge angles is observed except in the outer layer, where the changes in $Ma_{2}$, $Re_{2,\theta }$ and $T_w/T_2$ across the three cases may have an appreciable effect. Additional results from supersonic channel flow simulations by Modesti & Pirozzoli (Reference Modesti and Pirozzoli2016) at lower Mach numbers overlaid on figure 23 corroborate the common observation that the streamwise and wall-normal r.m.s. velocities close to the wall do not depend significantly on the Mach number when scaled with appropriate inner units.

Figure 29. (a) Streamwise and (b) wall-normal turbulence intensities at the streamwise location $(x-x_1)/{\delta _1^\star }= 580$ for $\alpha = 6^\circ$, $7^\circ$ and $8^\circ$, defined in terms of the density-weighted velocity scale $u_*=\sqrt {\overline {\tau _w}/\bar {\rho }}$, and plotted as a function of the semi-locally scaled wall-normal coordinate $y^\star = \bar {\rho } (\overline {\tau _w}/\bar {\rho })^{1/2}y/\bar \mu$. Also included is the reference data of channel case CH15C from Modesti & Pirozzoli (Reference Modesti and Pirozzoli2016) corresponding to $Re_{2,\tau } = 1015$ and $Ma = 1.5$.

Despite the reasonable agreements between WMLES and DNS outlined throughout this section, figure 30 shows that the core assumption of constant stress layer, represented by the momentum equation (B 1) of the equilibrium wall model, is not strictly satisfied within the wall-modelled layer for any of the three angles $\alpha =6^{\circ }$, $7^{\circ }$ and $8^{\circ }$. Specifically, the evaluation of the mean shear stress $\bar {\tau }$ provided in figure 30 using the DNS shows that the total stress $\bar {\tau }$ varies by amounts of order 10 % across the wall-modelled layer. These variations increase as the wedge angle decreases, or equivalently, as the friction Reynolds number decreases. In addition, as shown in figure 30, the ratio of the total and wall shear stresses, $\bar {\tau }/\overline {\tau _w}$, is not bounded by unity when the stresses are defined consistently with the Favre-averaged streamwise momentum equation. Instead, it features a maximum within the wall-modelled layer that was also observed in early computational work at lower Mach numbers by Gatski & Erlebacher (Reference Gatski and Erlebacher2002).

Figure 30. Direct numerical simulation wall-normal profiles of time- and spanwise-averaged shear stress $\bar {\tau }$ normalized with its value on the wall $\overline {\tau _w}$ at $(x-x_1)/\delta ^\star _1=580$ for wedge angles $\alpha =6^\circ$, $7^\circ$ and $8^\circ$. The dashed, dotted and solid lines denote the distributions of the viscous shear stress $\mu (\tilde {T}) (\partial \tilde {u}/\partial y)$, the turbulent shear stress $-\overline {\rho u''v''}$ and the total stress $-\overline {\rho u''v''} + \mu (\tilde {T}) (\partial \tilde {u}/\partial y)$. The vertical lines denote the matching location in the baseline WMLES (red lines) and in the coarse WMLES (blue lines).