2.1. Theoretical model

The three-dimensional unsteady Navier–Stokes equations are solved in strong conservation form. Upon non-dimensionalisation and mapping from Cartesian $(x,y,z)$ to generalised orthogonal curvilinear coordinates $(\xi ,\eta ,\zeta )$ , the governing equations are cast as

(2.1) \begin{equation} \frac {\partial }{\partial \tau }\Bigg (\frac {\boldsymbol{U}}{J}\Bigg ) + \frac {\partial {\boldsymbol{\hat {F}}}}{\partial {\xi }} + \frac {\partial {\boldsymbol{\hat {G}}}}{\partial {\eta }} + \frac {\partial {\boldsymbol{\hat {H}}}}{\partial {\zeta }} = \frac {1}{{\textit{Re}}_\infty }\Bigg ( \frac {\partial {\boldsymbol{\hat {F}}_v}}{\partial {\xi }} + \frac {\partial {\boldsymbol{\hat {G}}_{\!v}}}{\partial {\eta }} + \frac {\partial {\boldsymbol{\hat {H}}_{\!v}}}{\partial {\zeta }} \Bigg ) + \boldsymbol{S}, \end{equation}

where $\textbf {U} = [\rho \,\rho u \,\rho v \,\rho w\, E ]^T$ is the solution vector comprising conserved variables and $J$ connotes the determinant of the transformation Jacobian, given by $J = \text{det}(\partial {\xi _i}/\partial {x_j})$ . The terms $\boldsymbol{\hat {F}}$ , $\boldsymbol{\hat {G}}$ and $\boldsymbol{\hat {H}}$ signify the inviscid flux vectors in the $\xi ,\,\eta $ and $\zeta$ directions, respectively. Similarly, the viscous flux vectors in each coordinate direction are denoted by $\boldsymbol{\hat {F}}_v$ , $\boldsymbol{\hat {G}}_v$ and $\boldsymbol{\hat {H}}_v$ . For example, in the $\xi$ -direction, the flux vectors take the form

(2.2) \begin{gather} \boldsymbol{\hat {F}}=\left [\begin{array}{c} \rho \hat {U} \\ \rho u \hat {U}+\hat {\xi }_{x} p \\ \rho v \hat {U}+\hat {\xi }_{y} p \\ \rho w \hat {U}+\hat {\xi }_{z} p \\ ( E+p) \hat {U} \end{array}\right ] \,\,\,\,\,\,\,\,\text{and}\,\,\,\,\,\,\,\, \boldsymbol{\hat {F}}_{\!v} = \left [\begin{array}{c} 0 \\ \hat {\xi }_{x_i}\tau _{i1}\\ \hat {\xi }_{x_i}\tau _{i2}\\ \hat {\xi }_{x_i}\tau _{i3}\\ \hat {\xi }_{x_i}(u_j\tau _{\textit{ij}} - \varTheta _i) \end{array}\right]\!, \end{gather}

with the contravariant velocity $\hat {U} =\hat {\xi }_{x} u+\hat {\xi }_{y} v+\hat {\xi }_{z} w$ and the scaled metric terms $\hat {\xi }_i = J^{-1}\xi _i$ . Under the assumption of a perfect gas, the total energy variable may be derived from the primitive variables according to

(2.3) \begin{gather} E = \frac {p}{\gamma - 1} + \frac {1}{2}\rho (u^2 + v^2 + w^2). \end{gather}

The deviatoric stress tensor and the heat flux vector are given by

(2.4) \begin{gather} \tau _{\textit{ij}} = \mu \Bigg ( \frac {\partial {\hat {\xi }_k}}{\partial {x_j}}\frac {\partial {u_i}}{\partial {\hat {\xi }_k}} + \frac {\partial {\hat {\xi }_k}}{\partial {x_i}}\frac {\partial {u_j}}{\partial {\hat {\xi }_k}} - \frac {2}{3}\delta _{\textit{ij}}\frac {\partial {\hat {\xi }_l}}{\partial {x_k}}\frac {\partial {u_k}}{\partial {\hat {\xi }_l}} \Bigg ) \end{gather}

and

(2.5) \begin{gather} \varTheta _i = -\Bigg ( \frac {1}{(\gamma -1)M_\infty ^2}\Bigg )\frac {\mu }{\textit{Pr}_\infty }\frac {\partial {\hat {\xi }_j}}{\partial {x_i}}\frac {\partial {T}}{\partial {\hat {\xi }_j}} ,\end{gather}

respectively. In the above-mentioned formulations, both Stokes’ hypothesis for bulk viscosity $(\lambda = {-2}/{3}\mu )$ and Fourier’s law of heat conduction $(\varTheta _\xi \sim \partial T/\partial \xi )$ are invoked to provide closure for the viscous terms. The working fluid is assumed to be air, with Sutherland’s law used to impart temperature dependence to the viscosity coefficient and the Prandtl number held constant at a value of $\textit{Pr}_\infty =0.72$ . The flux vectors in the other computational coordinate directions are readily constructed in a similar fashion. Further theoretical details on casting the equations in the above-mentioned form may be found from Anderson et al. (Reference Anderson, Tannehill, Pletcher, Munipalli and Shankar2020). A vector of source terms, $\boldsymbol{S}$ , is retained in the governing equations to describe in a general manner any external excitation, injection, forcing or, as in the current effort, to enforce proper behaviour near the boundaries as discussed later.

2.2. Discretisation

Spatial discretisation is performed in a finite difference setting by constructing the flux vectors at each computational node and differentiating in the appropriate coordinate direction. For the inviscid fluxes in smooth regions of the flow and viscous fluxes throughout the domain, a non-dissipative sixth-order compact differencing scheme is employed to maintain high resolution across a broad range of spatial scales (Lele Reference Lele1992). An implicit eighth-order spatial filter is applied in each coordinate direction to provide regularisation to the solution, eliminating any spurious waves generated by parasitic dispersive behaviour at high wavenumbers (Visbal & Gaitonde Reference Visbal and Gaitonde2002). Following the recommendations of Schwartz (Reference Schwartz2023) for supersonic wall bounded flows, the free parameters of the filtering scheme are set to $\alpha = 0.45$ in the streamwise and wall-normal directions, and $\alpha = 0.49$ in the spanwise direction. In shocked regions of the flow, a third-order MUSCL type reconstruction is coupled with the approximate Riemann solver of Roe to compute the inviscid flux differences (Roe Reference Roe1981; van Leer & Nishikawa Reference van Leer and Nishikawa2021). Stability in the presence of shocks is further promoted by using the van Leer harmonic limiter, removing any Gibbs-type oscillations introduced by the polynomial reconstruction (LeVeque Reference LeVeque1992). As proposed by Visbal & Gaitonde (Reference Visbal and Gaitonde2005), successful hybridisation of the two schemes is accomplished by allowing an appropriate number of buffer points around shock waves for smooth transition from the compact differencing scheme to the upwind shock capturing scheme and back again. To effectively determine which regions of the solution require the use of shock capturing, a form of the modified Ducros sensor put forth by Larsson, Lele & Moin (Reference Larsson, Lele and Moin2007) is employed. The sensor is augmented with high-pass filter behaviour, using a formulation similar to that of Hendrickson, Kartha & Candler (Reference Hendrickson, Kartha and Candler2018) to prevent undesired dissipation of acoustic waves. High-order, nonlinearly stable time integration is achieved using the explicit, third-order SSP-RK3 algorithm (Shu & Osher Reference Shu and Osher1988).

2.3. Simulation configuration

The domain used for analysis is shown schematically in figure 2. Geometric dimensions are normalised using the channel half-height, $h$ , with the domain spanning $L_x \times L_y \times L_z = 28.5h \times 2h \times 2h$ . The streamwise length of the domain is fixed such that (i) prescribed synthetic inflow fluctuations discussed later recover to an equilibrium turbulent state before reaching the first shock wave; and (ii) the shock train generated by the enforced pressure rise achieves and then maintains an equilibrium position within the simulation domain. In the spanwise direction, sufficient extent is included for two-point correlations to decay. The domain is discretised into $1{,}639 \times 271 \times 195$ nodes in the streamwise, wall-normal and spanwise directions. For the first $25h$ , constant spacing is maintained in the streamwise direction. The final $50$ streamwise points are geometrically stretched by a factor of $1.03$ to provide an implicit sponge zone. Additional implicit damping by the mesh is included to improve the performance of the characteristic based boundary condition at the outlet (Colonius Reference Colonius2004).

Figure 2. Schematic of the computational domain. Dashed purple lines delineate spatial bounds of the explicit inflow sponge zone. The blue region indicates the spatial bounds of the implicit sponge zone where mesh stretching is applied.

Constant spacing is used in the spanwise directions. A hyperbolic tangent function is used in the wall-normal direction to both symmetrically distribute points about the centreline and cluster them near the walls. In the resolved portion of the domain, the streamwise spacing of the mesh is $\Delta x^+ = 14.1$ , based on the friction velocity, $u_\tau = \sqrt {\tau _w/\rho _w}$ , and viscous length scale, $\ell _v = \nu _w/u_\tau$ , at a reference station upstream of the interaction region $(x/h=8.5)$ . The spanwise spacing is $\Delta z^+ = 9.3$ , and the wall-normal spacing varies from $y^+_w = 0.85 \to y^+_\delta = 8.5 \to y^+_{cl} = 14.3$ at the wall, reference station boundary layer edge and duct centreline, respectively. Based on measured dissipation rates at the current flow conditions, the ratio of the effective mesh spacing $(\varDelta = \Delta x\boldsymbol{\cdot }\Delta y\boldsymbol{\cdot }\Delta z)^{1/3}$ to the local Kolmogorov length $(\eta = ({\nu ^3}/{\epsilon })^{1/4})$ does not exceed $8$ within the computed shock train.

The selected combination of discretisation schemes and mesh resolution allows the high-order filter to act as a sufficient energy drain and stabilise the calculation (Visbal & Gaitonde Reference Visbal and Gaitonde2002). Numerical damping is restricted to physical scales within the dissipative range of the wall-bounded flow (Jiménez Reference Jiménez2012) and modelling errors exist at only the smallest resolved and sub-grid scales. At present mesh resolutions, these errors have little impact on both solution accuracy (Kawai, Shankar & Lele Reference Kawai, Shankar and Lele2010; Garmann, Visbal & Orkwis Reference Garmann, Visbal and Orkwis2013) and measured turbulent transport properties (Morgan & Lele Reference Morgan and Lele2011). Consequently, no additional sub-grid scale model is invoked and the present computation is characterised as an implicit large eddy simulation.

To reduce the extent to which weak shock waves generated by the inflow boundary layer contaminate the interaction statistics, a sponge zone is used near the inflow plane to damp the conserved variables towards a set of reference values specified by the appropriate free stream conditions. Dashed purple lines in figure 2 denote the inflow sponge zone, which is situated in the region $(x,y) \in [0\,,2h]\times [({h}/{4})\,,({7h}/{4})]$ to allow for damping of undesirable shock waves. The location of the sponge is also chosen such that it does not extend into the inflow boundary layers, thereby leaving the prescribed inflow perturbations to develop properly. Attenuation strength decays as a function of both the streamwise and wall-normal coordinate, to provide a smooth transition between the damped solution inside the sponge region and the surrounding undamped solution. Further details on the class of inflow sponge used in this work are outlined by Morgan (Reference Morgan2012) and Huang (Reference Huang2022).

2.4. Initial conditions

The flow conditions for the simulation are patterned after the experimental campaign of Hunt & Gamba (Reference Hunt and Gamba2019), where a nominally Mach $2.0$ flow is supplied to a constant area, back-pressured channel. For computations, the Reynolds number based on the channel half-height is ${\textit{Re}}_h=4.3\times 10^5$ . This value is an order of magnitude lower than that of the experiments; as noted in several external (Priebe & Martín Reference Priebe and Martín2012; Adler & Gaitonde Reference Adler and Gaitonde2018) and internal flow investigations (Morgan et al. Reference Morgan, Duraisamy and Lele2014; Roussel Reference Roussel2016), this does not materially alter the scaled mean and instantaneous nature of the flow, but reduces mesh requirements and facilitates long integration times required for convergence of the measured turbulent statistics.

The flow is initialised by applying the Rankine–Hugoniot jump relations for a Mach $2.0$ normal shock at $L_x=3/4$ . Numerical stability is enhanced by this approach, as the initial shock waves formed in the domain are far removed from the computational outlet. After simulation start up, the normal shock in the initial condition naturally bifurcates and relaxes into the desired $\lambda$ -type shock train.

The simulation is marched from its initial conditions at a time step of $\Delta \tau = 0.001$ , with $\tau = t u_\infty /h$ being the non-dimensional time. In wall units, $\Delta \tau = 0.045\Delta t^+$ , which is well within the resolution range recommended for wall-bounded turbulence computations (Choi & Moin Reference Choi and Moin1994). Collection of flow statistics is initiated once the shock train attains a time invariant mean position within the domain, and is continued for $800$ non-dimensional time units $(h/u_\infty )$ , which translates to $32$ convective flow through times of the resolved portion of the domain and over $2200$ time units based on the boundary layer state upstream of the shock train.

2.5. Boundary conditions

Turbulent fluctuations are provided at the inflow plane via the synthetic digital filtering method described by Adler et al. (Reference Adler, Gonzalez, Stack and Gaitonde2018). First- and second-order moments attained from a precursor Reynolds-averaged Navier–Stokes (RANS) simulation are used to provide the required mean flow profiles and single point velocity correlations that generate the desired synthetic disturbances. Within the computational domain, the imposed perturbations interact with one another to produce turbulent boundary layers along the top and bottom walls which, after some development length, recover to an equilibrium turbulent state. To prevent spurious inflow pressure fluctuations from contaminating the internal flow interaction, the streamwise normal correlation is suppressed in the digital filtering routine (Ceci et al. Reference Ceci, Palumbo, Larsson and Pirozzoli2022; Laguarda & Hickel Reference Laguarda and Hickel2024).

To maintain proper boundary layer behaviour throughout the duct, the no-slip condition is enforced on the three velocity components along the top and bottom walls. The walls are also treated as iso-thermal, with the wall temperature given by the recovery temperature of the free stream flow, $T_r = T_\infty (1+ r({(\gamma - 1)}/{2})M_\infty ^2)$ , and $r = \textit{Pr}_\infty ^{1/3}$ . The flow is assumed homogeneous in the spanwise direction, allowing use of periodic boundary conditions. While the use of periodicity makes a direct comparison with the reference experimental results difficult, as shown in both previous shock train computational studies (Morgan et al. Reference Morgan, Duraisamy and Lele2014; Roussel Reference Roussel2016) and the present work (§ 3.4), the key physics of shock trains are still captured under these simplified conditions.

Along the outflow plane, an elevated back-pressure at $3.4$ times the mean inflow pressure (i.e, $\overline {p}_b=3.4\overline {p}_i$ ) is imposed using a nonlinear, partially non-reflecting characteristic boundary condition (Thompson Reference Thompson1990; Poinsot & Lele Reference Poinsot and Lele1992). The characteristic formulation allows flow disturbances to exit the computational domain while maintaining the outlet pressure near the desired value. To prevent non-physical coupling between incoming acoustic modes and outgoing vorticity and entropy modes, the partially non-reflecting condition is augmented with multi-dimensional treatment (Yoo & Im 2007; Lodato, Domingo & Vervisch Reference Lodato, Domingo and Vervisch2008; Fosso et al. Reference Fosso, Deniau, Lamarque and Poinsot2012), which minimises spurious acoustic reflection from the outlet in the presence of the turbulent, subsonic flow expected downstream of the shock train interaction. The damping provided by the geometric stretching further improves the behaviour of the computational outlet (Colonius Reference Colonius2004). To more accurately replicate the physics occurring within a shock train, the partially non-reflecting condition is applied only to the subsonic portions of the boundary, i.e. the local acoustic field is modified only where the local flow is subsonic. In locally supersonic regions, the evolution of the acoustic field is determined entirely by the flow within the computational boundary.

The pressure histories of the inflow and outflow boundaries are investigated to confirm that the boundary condition formulations do not introduce spurious frequency content. Figure 3 contains time and frequency domain descriptions of the wall pressure at the inflow and outflow boundaries. In the time domain (figure 3 a), it is clear that the desired inflow and outflow pressures are being enforced. Power spectral densities, $\phi (St)$ , of the outlet pressure (figure 3 b), of the inflow and outflow pressure signals demonstrate no spurious low frequency, $St \leqslant \mathcal{O}(.01)$ , that might couple with the inherent low-frequency motion observed in shock trains (Hunt & Gamba Reference Hunt and Gamba2019). The slight peak observed at $St\sim 0.1$ in the outlet pressure spectrum is due a naturally occurring undulation of the shock train and corresponds to a minor mode of oscillation that has also been observed experimentally (Hunt & Gamba Reference Hunt and Gamba2019).

Figure 3. (a) Time domain and (b) frequency domain descriptions of the inflow and outflow boundaries. No spurious spectral peaks, $St\leqslant \mathcal{O}(.01)$ , that might drive the actual interaction are observed.

3.1. Instantaneous flow structure

The salient features of the interaction are highlighted in figure 4 by considering an instantaneous snapshot of the full three-dimensional shock train flow. Dilatation isosurfaces, coloured by local static pressure, identify the series of shock waves that emerge within the channel. The leading shock is bifurcated, displays a $\lambda$ structure and retains a short Mach stem. Subsequent shock waves do not bifurcate and retain an orientation approximately normal to the local flow. To visualise the instantaneous boundary layer structure, an isosurface of Q-criterion (Jeong & Hussain Reference Jeong and Hussain1995) coloured by streamwise velocity is employed. In the upstream region of the shock system ( $x/h\leqslant 12$ ), the wall bounded turbulent flow is densely populated by arrays of vortex filaments, which are particularly visible in the local outer layer. At the foot of the leading STBLI ( $x/h\sim 14$ ), the boundary layer decelerates and is deflected away from the wall, resulting in a thicker near-wall turbulent layer that persists for the remaining streamwise extent of the domain. Subsequent interactions between normal shock waves and the near-wall layer are readily identified. Mach number contours show that the shock train induces a departure of the sonic line away from the wall, and the flow attains a mixed supersonic and subsonic turbulent state at the end of the resolved portion of the computational domain.

Figure 4. Instantaneous visualisation of the full 3-D, spanwise periodic flow field. Shown are an isosurface of Q-criterion coloured by streamwise velocity (near-wall vortical structures), an isosurface of dilatation coloured by static pressure (grey geometric surfaces), and contours of the Mach field ( $x{-}y$ plane at the left channel face and $y{-}z$ plane at the channel outflow).

3.2. Interaction features and data collection

To facilitate description of the evolution of the flow statistics, both globally across the entire shock train and locally due to each individual STBLI shown in figure 4, it is convenient to formally parse the shock train and split it into various ‘regions’. The mean wall pressure is a suitable variable for this purpose, with each STBLI region bounded on either side by a minimum in the streamwise wall pressure gradient, shown in figure 5. Delineating the train into individual STBLIs allows for an analysis of each interaction with a satisfactory degree of isolation and aids in determining how global flow changes across the shock train emerge as the aggregate of localised changes.

Upstream of the shock waves, the turbulent boundary layers (TBLs) display self-similar, equilibrium turbulent profiles and are free from the upstream influence of the leading STBLI. This area, referred to as the ‘upstream region’, extends from $8 \leqslant x/h \leqslant 12$ . The region where the first bifurcated STBLI resides ranges from $12 \lt x/h \leqslant 15.8$ , and is referred to as ${R}_1$ . The zone prescribed to the second STBLI is located in the range $15.8 \lt x/h \leqslant 17.5$ , denoted as ${R}_2$ . The regions of the third $({R}_3)$ and fourth $({R}_4)$ STBLIs range over $17.5 \lt x/h \leqslant 18.9$ and $18.9 \lt x/h \leqslant 20$ , respectively. Finally, the area between the upper bound of ${R}_4$ and the end of the resolved computational domain $(x/h=25)$ is referred to as the ‘downstream’ region. For ease of reference, a tabulated summary of the defined region bounds is given in table 1.

Figure 5. Mean streamwise wall pressure and its gradient. Open symbols denote the local minima in wall pressure gradient, used to define spatial bounds on each individual STBLI. Vertical dotted lines are included to readily delineate STBLI region boundaries.

Table 1. Description of described shock train regions.

While more than four distinct STBLIs are apparent in the channel (figure 4), it is sufficient for clarity and brevity to examine details of the first four interactions. Beyond ${R}_4$ , as shown in subsequent sections, the local Mach number approaches unity and shocks have a comparatively weak effect on the mean flow statistics.

Several data collection procedures are used to accurately assess the modification of the mean flow structure and near-wall turbulence by the overall shock train and individual STBLIs. Rolling averages of the primitive variables, velocity covariances and density weighted velocity covariances are collected at each grid point over a time interval of $800\tau$ . Time-resolved velocity signals are collected across two spanwise lines at fixed height to compute the evolution of correlation coefficients. At the interfaces between shock train regions ( $x/h\in \{12,\,15.8,\,17.5,\,18.9,\,20,\,25\}$ ), $y{-}z$ planes are extracted at a sampling rate of $St = fh/u_\infty =25$ to enable resolved computation of turbulent kinetic energy budgets and energy wavenumber spectra. In addition to being averaged in time, flow statistics are averaged in the spanwise direction wherever appropriate, to significantly increase the available sample size. Note that throughout the analysis of the turbulent interaction, both Reynolds $(\phi = \overline {\phi } + \phi ')$ and Favre $(\phi = \widetilde {\phi } + \phi^{\prime\prime} ,\,\widetilde {\phi } = \overline {\rho \phi }/\overline {\rho } )$ decompositions are used to assess flow behaviour.

3.3. Characterisation of the incoming turbulent flow

To assess the resolution of the computational mesh and validate the nature of the incoming turbulent boundary layers, velocity statistics of the flow in the upstream region are examined in detail. Data are extracted at the streamwise station $x/h=8.25$ , as this location provides enough development length for the synthetic inflow perturbations to recover to an equilibrium turbulent state (Wu Reference Wu2017) and possesses local Reynolds numbers $({\textit{Re}}_\tau ,{\textit{Re}}_\theta ,{\textit{Re}}_{\delta _2})$ that compare favourably with reference DNS data used to validate the turbulent statistics. Table 2 lists key properties, including edge Mach number, boundary layer thicknesses, relevant Reynolds numbers based on boundary layer properties, skin friction coefficient, incompressible shape factor and peak turbulent Mach number.

Table 2. Local turbulent boundary layer properties at the reference station $x/h=8.25$ . ${\textit{Re}}_\tau = \delta /\ell _v,\,{} {\textit{Re}}_\theta = \rho _\infty u_\infty \theta /\mu _\infty ,\, {\textit{Re}}_{\delta _2} = \rho _\infty u_\infty \theta /{\mu _w},\, C_{\!f} = 2\tau _w/\rho _\infty u_\infty ^2, \,H_i = \delta ^*_i/\theta _i,\,M_t=\sqrt {\overline {u_i'u_i'}}/\overline {c}$ .

Given the difficulties with defining a boundary layer ‘edge’ location in the classical sense for the present case comprising internal flow STBLIs, the local boundary layer edge ( $\delta )$ is defined using a reference vorticity method (Pasquariello, Hickel & Adams Reference Pasquariello, Hickel and Adams2017). In this approach, $\delta$ is measured using its classical definition at some undisturbed upstream station. Vorticity at the boundary layer edge is then assumed to be constant over all downstream boundary layer edge locations and used as a threshold to identify the wall normal location at which the TBL terminates. At locations upstream of the shock train interaction, this definition relaxes to the traditional one of boundary layer thickness while providing a more robust description of the layer edge location in the presence of STBLIs. Standard definitions for compressible and incompressible integral thicknesses are used to calculate boundary layer properties.

Figure 6. Spanwise correlation coefficient of the three velocity components over the half-span of the computational domain ( $0 \leqslant z/h \leqslant 1$ ), at the reference station $x/h=8.25$ . Open symbols correspond to data recorded at location of peak streamwise velocity variance ( $y^+ \sim 11$ ). Closed symbols correspond to data recorded at local boundary layer edge ( $y=\delta$ ).

To determine if the domain retains an appropriate spanwise length, the correlation coefficient of the three velocity components is examined. Figure 6 shows the correlation curves at two wall-normal heights, one at the inner layer streamwise velocity variance peak ( $y^+\sim 11$ ) and one at the local boundary layer edge ( $y = \delta$ ). Inner layer spanwise correlations decay to negligible values within approximately $100\ell _v$ , coinciding with the approximate viscous spacing of near-wall streaks. Outer layer correlations decay over a larger spatial scale, dropping to negligible values over lengths of $\mathcal{O}(\delta )$ . As no unnatural statistical correlations are observed over the spanwise direction, the spatial extent of the domain is deemed to be of sufficient extent to preclude imposition of spurious wavelengths.

The integral transform of Van Driest,

(3.1) \begin{equation} u_{vd}=\int ^{u^+}_0 \bigg (\frac {\overline {\rho }}{\overline {\rho }_w}\bigg )^{(1/2)} \,{\rm d}u^+ ,\end{equation}

is also applied to the mean streamwise velocity profile at the reference station. The resulting distribution is reported in figure 7, where it is plotted against the compressible Mach 2 DNS data of Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011) ( ${\textit{Re}}_\tau =250,{\textit{Re}}_{\delta _2} = 715,{\textit{Re}}_\theta = 1122$ ) and the incompressible DNS data of Schlatter & Örlü (Reference Schlatter and Örlü2010) ( ${\textit{Re}}_\tau = 270,{\textit{Re}}_{\delta _2} = 670$ ). For visualisation, only half of the wall-normal profile from the current calculation has been considered, as the time-averaged variations of velocity and density are symmetric about the channel centreline. Good agreement is observed between the present transformed compressible profile and the reference equivalent incompressible velocity profiles. When presented in inner unit scaling, each profile exhibits excellent collapse in the near-wall linear sublayer. At the considered Reynolds number, approximately logarithmic scaling is present for a brief interval inside the overlap region in the range $30 \leqslant y^+ \leqslant 100$ . As $y^+$ continues to increase, the profiles adopt a wake-like structure, representative of the outer boundary layer flow.

Figure 7. Mean streamwise velocity profile plotted in equivalent inner units, extracted from the reference location at $x/h=8.25$ . Open circles correspond to the Van Driest transformed data of Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011). Closed triangles correspond to the inner unit data of Schlatter & Örlü (Reference Schlatter and Örlü2010). Here, $k=0.41$ and $c = 5.1$ are used for the reference log-layer scaling of $u^+(y^+) = ({1}/{k})\text{log}(y^+) + c$ .

Wall-normal profiles of the single point velocity correlations also compare favourably to those from compressible DNS data. Alignment between the present simulation and the reference data is maintained when examining the flow in both inner (figure 8 a) and outer (figure 8 b) scaling. In particular, the inner scaling highlights that the peak of streamwise velocity fluctuations occurs at approximately $y^+=11$ , consistent with the well-documented physics of wall bounded turbulent flows (Smits, McKeon & Marusic Reference Smits, McKeon and Marusic2011). The slight disagreement between the velocity fluctuations in the current calculation and the comparison DNS are not consequential and may be attributed to small variations in Reynolds number and mesh resolution between the two databases (Poggie, Bisek & Gosse Reference Poggie, Bisek and Gosse2015).

Figure 8. (a) Inner unit scaling. (b) Outer unit scaling. Wall-normal profiles of single point velocity correlations at the reference station $x/h=8.25$ , shown in both inner and outer scaling. Line plots denote present LES data, open symbols denote reference DNS data of Pirozzoli & Bernardini (Reference Pirozzoli and Bernardini2011). Profiles have been normalised by the square of the friction velocity at the boundary layer reference station.

Figure 9. Comparison of the wall pressure profiles in the current LES and the reference experiments of Hunt & Gamba (Reference Hunt and Gamba2019). Note the profiles have been shifted in the streamwise direction from their original locations such that their initial pressure rises are coincident.

3.4. Validation of shock train behaviour

In addition to characterising the state of the incoming turbulence, it is necessary to examine whether or not the present simulation faithfully reproduces the expected dynamics of the shock train interaction. Time-averaged wall pressure distributions are compared between the present computation and the experimental data of Hunt & Gamba (Reference Hunt and Gamba2019) in figure 9. Favourable agreement is apparent between the two normalised profiles both upstream of and during the initial wall pressure rise. Modest deviation is apparent in the more downstream portion of the shock train, with the numerical pressure rise occurring more gradually than the experimental one. Given that the mean wall pressure is both a measure of the compressive processes occurring due to outer layer shock waves and the non-local contributions from boundary layer turbulence (Smits & Dussauge Reference Smits and Dussauge2006), the distributions in figure 9 indicate that in a mean sense the critical physics of the experimental shock train are being replicated numerically.

Figure 10. Comparison of the normalised spectra of streamwise shock wave oscillations between the present LES and the reference experimental data set. (a) Leading leg of first bifurcated shock. (b) Trailing leg of first bifurcated shock. (c) Second shock in train. (d) Third shock in train. In all figures, the frequency axis is scaled using the characteristic frequency of the channel half-height ( $St_h$ ). Consistent with experiments, normalisation has been carried out using the variance of shock position fluctuations ( $\sigma ^2(x'_i)$ ) and the fundamental frequency of the shock train $f_c=L_{st}/U_\infty$ , where $L_{st}$ is the characteristic length of the full shock system.

To ensure dynamic processes are also being captured by the simulation, statistics related to shock wave oscillations are presented. For this purpose, we examine the spectral behaviour of shock wave motions in the present LES, as streamwise oscillations of the shock system are a characteristic feature of shock train interactions (Matsuo et al. Reference Matsuo, Miyazato and Kim1999; Gnani et al. Reference Gnani, Zare-Behtash and Kontis2016). The power spectral densities (PSD) of four shock features are presented in figure 10. The wave features being considered are the leading leg of the first bifurcated shock (figure 10 a), the trailing leg of the first bifurcated shock (figure 10 b), the shock of the second STBLI in the train (figure 10 c) and the shock of the third STBLI in the train (figure 10 d). These shocks are labelled $x_1,x_2,x_3$ and $x_4$ respectively. Oscillatory behaviours of each shock are also compared with the equivalent wave features taken from the experiments of Hunt & Gamba (Reference Hunt and Gamba2019).

Figure 11. Streamwise pressure profiles from the shock train interaction. Curves are the centreline ( $p_{cl}$ ) and wall ( $p_w$ ) profiles from the LES. Open symbols denote pressure values computed using Fanno flow relations.

All shock features display similar spectral slopes to their experimental counterparts, characterised by broadband oscillations with high energy encoded at larger time scales and an approximate power law decay as frequency increases. Furthermore, minor modes of shock oscillation near $St_h\sim 0.1$ that appear in experiments are accurately captured in the computational spectra. These are most noticeable in downstream shock waves (see shock $x_4$ , figure 10 d), where harmonic, undulating motions of the shock system become more pronounced. Agreement between computation and experiment in both the dynamic motion of the strongest shock waves in the system and the time mean compression is encouraging, and indicates that, even under the simplified conditions of spanwise periodicity, proper shock train behaviour is observed in the present calculation.

4.1. Mean interaction structure

Streamwise pressure distributions, which exemplify competing mechanisms of compression in supersonic and subsonic regions of the shock train, are plotted both along the bottom wall and duct centreline in figure 11. The wall and centreline pressures profiles collapse upstream of the interaction, indicating negligible wall normal pressure gradients and an essentially one-dimensional behaviour of the flow upstream of the shock train. A weak adverse pressure gradient is measured in this region (figures 5 and 11). Given the nature of the flow upstream of the shock train, it is expected that a quasi-one-dimensional model might predict the initial compression occurring in the channel. To corroborate this idea, Fanno flow theory is invoked (Shapiro Reference Shapiro1983). Using a friction factor based on the average upstream skin friction coefficient $(\bar {f} = 0.04)$ , the pressure profiles in the channel configuration before the interaction are well predicted by classical theory, connoting that only near-wall boundary layer processes alter the mean pressure distribution upstream of the shock train.

Once the flow reaches the shock train, pressure distributions between the walls and the centreline begin to diverge. The wall pressure exhibits a smooth and monotonic rise from the pressure directly upstream of the shock train to the prescribed time-averaged outflow pressure. In contrast, sharp changes emerge in the streamwise pressure profile along the centreline as a result of rapid compression via shock waves. Steep increases in pressure along the centreline are followed by subsequent decreases due to reacceleration of the flow, which are consequences of streamline curvature caused by boundary layer growth and expansion fans emanating from the locations where shock waves impinge on the sonic line near the walls (Howarth Reference Howarth1948; Geerts & Yu Reference Geerts and Yu2016).

After the first STBLI, the near-wall subsonic flow has not been compressed to the prescribed downstream pressure. The process of flow deceleration, boundary layer thickening and streamline deflection occurs again downstream of the first shock structure, leading to the formation of a second STBLI. This process repeats, until many STBLIs are formed in the channel and the subsonic flow reaches the prescribed downstream pressure.

To study the decelerating nature of the shock train, the mean Mach number field is shown in figure 12 in the form of two-dimensional contours. The bifurcated nature of the leading shock is clearly evident, with a leading leg characterised by a fan of oblique compression waves, typical of impinging interactions, while the trailing leg is nearly normal to the local flow and induces rapid compression (Green Reference Green1970). Near the channel centreline, a short Mach stem is also present in the time-averaged structure of the leading shock. Subsequent downstream shock waves maintain a near-normal orientation to the flow, also engendering intense compression, demarcated by sharp transitions from supersonic to subsonic values. The sonic line, shown in black within the figure, initially remains close to the wall upstream of the shock train before detaching once the first STBLI is reached. As a result, a large portion of the flow subjected to adverse pressure gradients is subsonic, displaying a more gradual and fully continuous compression until the applied back pressure is attained.

Figure 12. Time- and spanwise-averaged Mach number field, $M(x,y)$ , for the shock train interaction. Black lines indicate the time-mean $M=1$ isosurface, dividing subsonic and supersonic regions of flow.

Further detail on the local variation of Mach number is examined through its variation in the streamwise direction. Figure 13 shows $M(x)$ curves at three wall-normal locations; $y=h$ , $y=h/2$ and $y=h/4$ . The centreline profile, $y=h$ , demonstrates rapid inviscid compression and expansion across each STBLI. At further downstream locations, the centreline flow is compressed by intermittent and weak shocks, resulting in smoother time mean curves. Past the strongest interactions in the train ( $R4$ and beyond), the centreline Mach number remains above unity for the remainder of the interaction. Channel quarter height and eighth height profiles undergo more gradual compression through the leading STBLI, and transition from supersonic to subsonic values as the boundary layer thickens across the leading oblique leg. Note these profiles occur at $y=h/2$ and $y=h/4$ , respectively, as the full channel height is defined as $L_y=2h$ . Oscillations of these three profiles decay downstream, indicating that as the core Mach number decreases, the acceleration of the outer inviscid flow exerts a weaker influence on the underlying, perpetually subsonic flow.

Figure 13. Streamwise distributions of Mach number at three wall-normal locations, $y=h$ , $y=h/2$ and $y=h/4$ .

The structure of the streamwise velocity field further demonstrates the undulating effect induced by repeated inviscid acceleration and deceleration. Isovels in figure 14 depart from the wall and remain detached after the first STBLI, with associated fluctuation of the boundary layer edge and sonic line heights, which are overlaid on the contours. The first STBLI notably increases the boundary layer thickness from its approach value which is approximately $\delta /h\sim 0.4$ . As discussed in more detail in § 5, the line of maximum TKE is also deflected away from the wall. Subsequent STBLIs locally modulate the boundary layer edge, with shock waves thickening the layer and expansion fans thinning the layer.

Figure 14. Mean streamwise velocity field. Note velocity has been normalised by the free stream velocity, thus the colour bar indicates isovels of $u/u_\infty$ . The dashed white line indicates the local boundary layer edge. The dashed black line denotes the locus of points where turbulent kinetic energy attains a local maxima, which is used as an approximate identifier for the midpoint of the thickened wall layer.

The coupled effect of the streamline curvature also introduces strong oscillating vertical velocity components in the flowfield. Contours of $v/u_\infty$ in figure 15 demonstrate the prominence of these $\pm v$ variations in the local outer boundary layer, and also highlight the observation that the wall-normal velocities decrease in magnitude as the local shock waves decrease in strength. The variations in vertical velocity introduce new outer layer mass and momentum conservation effects into the turbulent flow that are absent in the upstream, undisturbed boundary layer. Furthermore, the imposition of vertical velocities on the outer edge of the flow may be thought of as local $v(x)$ type boundary conditions on the underlying turbulence layer, facilitating analogies to be drawn between the present complex shock train case and the more canonical problem of an incompressible boundary layer separating due to a specified suction/blowing function above the wall (Kitsios et al. Reference Kitsios, Sekimoto, Atkinson, Sillero, Borrell, Gungor, Jiménez and Soria2017; Wu, Meneveau & Mittal Reference Wu, Meneveau and Mittal2020).

Figure 15. Mean vertical velocity field, $v/u_\infty$ . Overlaid curves retain their definitions from figure 14 .

Figure 16. Streamwise profiles of boundary layer and displacement thicknesses. Profiles have been normalised by the value at the location upstream boundary of ${R}_1$ , which is free from the upstream influence of the leading STBLI. Open circles denote the thicknesses predicted by the ${\textit{Re}}_x^{-1/7}$ power law of Prandtl.

Streamwise distributions of the boundary layer and displacement thickness quantitatively describe this trend. Figure 16 shows these two thicknesses, normalised by the value at the leading edge of ${R}_1$ . Upstream of the shock train, the boundary layer thicknesses are predicted by Prandtl’s $1/7$ power law, given by

(4.1) \begin{equation} \delta \sim x {\textit{Re}}_x^{-1/7}. \end{equation}

As each interaction is crossed, the boundary layer and displacement thickness undergo a net increase, attaining a local maximum before decreasing as the end of an STBLI region is approached. Notably, the displacement thickness exhibits a near four-fold increase under the first STBLI, then fluctuates at decreasing amplitudes before asymptotically approaching a value slightly lower than $3.5$ times its value at the reference station. Significant variations in the displacement layer thickness are consistent with the ‘aerodynamic nozzle’ hypothesis of spatially confined STBLI flows (Geerts & Yu Reference Geerts and Yu2016), where streamline deflection, induced by near-wall aerodynamic blockages, is in part responsible for the re-acceleration of the flow in the inviscid core of the channel.

4.2. Reverse flow structure

Depending on the strength of individual STBLIs in the system, significant local changes in the mean wall shear are expected. To evaluate these effects, the Van Driest II transformation (Hopkins & Inouye Reference Hopkins and Inouye1971) is used:

(4.2) \begin{equation} C_{f_i} = F_cC_{\!f},\,\,\,\,\,\,{\textit{Re}}_{\theta _i}=F_\theta {\textit{Re}}_\theta, \end{equation}

wherein

(4.3) \begin{equation} F_c = \frac {\overline {T}_w/T_\infty - 1}{\arcsin ^2\alpha },\quad F_\theta = \frac {\mu _\infty }{\overline {\mu }_w} ,\quad \alpha = \frac {\overline {T}_w/T_\infty - 1}{\sqrt {\overline {T}_w/T_\infty (\overline {T}_w/T_\infty - 1)}}. \end{equation}

The incompressible correlation

(4.4) \begin{equation} C_{f_i} = 0.026/{\textit{Re}}_{\theta _i}^{1/4} \end{equation}

is overlaid on the transformed skin friction curve in figure 17. The results show that the skin friction relaxes to equilibrium values upstream of the interaction. The distribution of computed skin friction in figure 17 also clearly shows that wall shear is altered by individual STBLIs, with skin friction attaining local minima in regions $R1\to R4$ . Further downstream, minima are readily attributed to the weaker STBLIs in the channel.

Figure 17. Streamwise distributions of equivalent incompressible skin friction coefficient. Solid line denotes the computed mean shear, dashed line corresponds to the Blasius skin friction correlation of (4.4).

The leading STBLI is strong in the sense that it induces time mean separation, while successive downstream interactions do not. Since the local outer Mach number decreases with streamwise distance, the successively weaker shocks have correspondingly milder influence on the local skin friction. Soon after the end of ${R}_4$ , the impact of the mean flow shock waves on the wall skin friction becomes practically negligible.

A probabilistic description of the reverse flow is presented in figure 18, by evaluating for all $x/h$ the percentage of time, $\gamma (x)$ , backflow is observed just off the wall. To describe in detail the statistical nature of the flow separation, we briefly adopt the nomenclature of Simpson (Reference Simpson1989), who, in the context of two-dimensional turbulent separations, described three classes of backflow regions. These are transitory detachment (TD), intermittent transitory detachment (ITD) and incipient detachment (ID), which respectively correspond to reverse flow probabilities of $\gamma \geqslant 50\,\%$ , $50 \,\%\,\gt \,\gamma \geqslant 20\,\%$ and $\gamma \leqslant 1\,\%$ . Under this framework, the first STBLI in the shock train induces two regions of transitory detachment. These regions correspond to the mean separation observed in figure 17. The next three interactions produce intermittent transitory detachment (corresponding to incipient separation), with back flow probabilities in excess of 20 %. Further downstream, interactions are too weak to meet this requirement and the reverse flow probability gradually approaches incipient detachment as the streamwise coordinate increases. Returning to the terminology of STBLIs, the shock train under investigation can be classified as having a strong leading STBLI (inducing mean separation), followed by three incipient interactions, which in turn are followed by even weaker interactions.

Figure 18. Likelihood (as a percentage) of negative shear stress along the bottom wall. TD, transitory detachment; ITD, intermittent transitory detachment; ID, incipient detachment.

Figure 19. Enlarged view of turbulent separation bubble structure, with the wall-normal coordinate normalised by the upstream reference viscous length, $\ell _{v,{\boldsymbol{ref}}}$ . The chained black line denotes the location of local maxima in TKE. The solid white line denotes the boundary of reverse flow regions, i.e. the set of off-wall points such that $\overline {u}=0$ .

The two-dimensional structure of the train is used to describe the collapsed nature of the leading mean separation bubble. Contours of streamwise velocity are plotted in figure 19 on a highly distorted axis scale for illustration. Two distinct regions of mean back flow are observed, with the first region maintaining a thin height. When scaled by the upstream viscous length, the height of the leading bubble is calculated to be $ y^+_{\textit{ref}}\sim 5$ . Viscous bubble structures such as the one observed are characteristic of low Reynolds number interactions, as noted by Larsson et al. (Reference Larsson, Kumar, Oberoi, Renzo and Pirozzoli2022) and Laguarda et al. (Reference Laguarda, Hickel, Schrijer and Van Oudheusden2024). A quantitative description of the bubble separation and attachment locations is given in table 3.

Table 3. Time mean separation (S) and reattachement (R) locations.

4.3. Emergence of self similarity

The presence of rapid undulations in the pressure field breaks the notion of mechanical equilibrium and self-similarity in wall-bounded flows (Clauser Reference Clauser1954). However, in the limit of these pressure variations tending to very small values, self-similar states can manifest even under conditions of non-zero pressure gradient (Devenport & Lowe Reference Devenport and Lowe2022). Evidence that this behaviour can persist in multi-shock systems was recently presented by Zuo et al. (Reference Zuo, Memmolo, Huang and Pirozzoli2019) for the case of conical shock generator/flat plate supersonic interaction. In the present configuration, jump discontinuities become progressively weaker across each individual shock in the train, and the variations of the mean streamwise pressure gradient asymptotically decay (figure 5), potentially allowing a self-similar state to arise at far-downstream locations.

Figure 20. Streamwise velocity profiles in outer coordinates with global scaling. Open circles on each profile indicate the computed location of the boundary layer edge.

To determine if this is the case for the present problem, the evolution and scaling of the mean velocity profile in the shock train is investigated. Velocity profiles, presented in a ‘global’ outer unit scaling, are shown in figure 20. Here, the profiles are normalised by the defining length and velocity scales of the entire flow, which are respectively the channel half-height $(h)$ and the incoming free stream velocity $(u_\infty )$ . Upstream of the interaction, at the recorded $x/h=12$ and $x/h=13.45$ locations, the streamwise velocity profile is full, consistent with the observed equilibrium boundary layer behaviour observed ahead of the interaction. At all subsequent downstream stations, the mean streamwise velocity profiles maintain an inflection point. A qualitatively similar inflectional profile was noted for the transonic diffuser interactions explored by Sajben et al. (Reference Sajben, Morris, Bogar and Kroutil1991), where flow in a shocked transonic diffuser was subsequently exposed to a mild downstream adverse pressure gradient. The extended presence of a mild adverse pressure gradient in the downstream region of the current shock train induces similar pressure gradients on the near-wall flow, indicating that such a velocity profile is common for these types of flow conditions.

To better assess if self-similarity emerges as the flow recovers downstream, a change to ‘local’ boundary layer coordinates is invoked. As noted by Castillo & George (Reference Castillo and George2001), use of an appropriate outer scaling should consolidate velocity profiles that are in a self-similar adverse pressure gradient state. Such a scaling is given by the defect law of Zagarola & Smits (Reference Zagarola and Smits1998), who took the outer velocity $(u_o)$ scale to be $u_o=u_e ({\delta ^*}/{\delta })$ and the length scale to be $y/\delta$ . Distributions resulting from this variable transform are shown in figure 21.

Figure 21. Streamwise velocity profiles in the local defect scaling. Under this scaling, the deviation from a canonical full profile and adoption of an inflected self-similar profile becomes apparent.

Profiles from upstream of the interaction exhibit appreciable collapse under this scaling, though the upstream influence of the leading STBLI results in some deviations between the two profiles before the separation bubble. In the presence of the strong leading STBLI, poor agreement is observed. However, near subsequent interactions, extended regions of outer layer agreement are observed. Further downstream, profiles exhibit much better collapse, indicating the turbulent flow is converging to the self-similar profile that might be expected of an adverse pressure gradient equilibrium state. Gradual relaxation to this state is expected to be a consequence of the diminishing pressure gradient imposed by the unique outer shock structure of the shock train interaction, as seen in the recovery behaviour of turbulent boundary layers in single transonic interactions (Sajben et al. Reference Sajben, Morris, Bogar and Kroutil1991; Pirozzoli et al. Reference Pirozzoli, Bernardini and Grasso2010).

Figure 22. Equilibrium evolution of the observed turbulent boundary layer, plotted using the self-similarity coordinates of Castillo & George (Reference Castillo and George2001). Dashed lines denote zero and adverse pressure gradient relationships between the two coordinates. The inset contains chained lines to help visualise the path that the flow takes through the variable space. White triangles indicate bounds of ${R}_1\to {R}_4$ , with the final white triangle denoting the final streamwise station of the resolved portion of the domain.

The departure and return to self-similarity in the velocity profiles indicates that the boundary layer may be transitioning from one equilibrium state to another. The equilibrium parameter of Castillo & George (Reference Castillo and George2001) is used to examine this eventuality. Upon introducing the pressure gradient parameter in outer layer variables,

(4.5) \begin{equation} \varLambda = \frac {\delta }{\rho _eu_e^2{\rm d}\delta /{\rm d}x}\frac {{\rm d}p_w}{{\rm d}x} = -\frac {\delta }{u_e {\rm d}\delta /{\rm d}x}\frac {{\rm d}u_e}{{\rm d}x}, \end{equation}

the authors suggest that most wall boundary layers are in an equilibrium state, under the condition $\varLambda = \mathrm{const}$ . in the flow. Within this constraint, profiles of streamwise velocity may be collapsed since

(4.6) \begin{equation} \delta \sim u_e^{-1/\varLambda }. \end{equation}

Three distinct values of $\varLambda$ exist, one each for the case of favourable, zero and adverse pressure gradient boundary layers (FPG, ZPG, APG), with $[\varLambda _{\textit{FPG}},\,\varLambda _{\textit{ZPG}},\,\varLambda _{\textit{APG}}] = [-1.92,\,0,\,0.22]$ . For the present case, by plotting $u_e/u_\infty$ versus $\delta /\delta _{in}$ in figure 22, the existence (or lack thereof) of equilibrium turbulent states may be verified. In the upstream portion of the interaction, the flow is approximately in a zero pressure gradient equilibrium state, but begins to depart as the flow approaches the first interaction. Between $R1$ and $R4$ , the turbulent flow is substantially out of equilibrium, but as the outlet of the domain is approached, the flow asymptotically approaches the predicted adverse pressure gradient equilibrium. The inset of figure 22 contains a connected line plot demonstrating the path travelled in $u_e/u_\infty -\delta /\delta _{in}$ space, highlighting how shocks induce large variations in boundary layer edge location and velocity within each STBLI region. However, these deviations diminish in strength as the end of the computational domain is approached, and the flow gradually relaxes to the predicted APG equilibrium trajectory.

The unique outer layer shock structure of the internal flow shock train interaction results in the state of the near-wall boundary layer taking a starkly different trajectory through the $u_e/u_\infty -\delta /\delta _{in}$ space compared with other external flow interactions. In supersonic impinging cases with a single shock, the relation between boundary layer growth and edge velocity relaxes to the equilibrium upstream zero pressure gradient condition (Yu et al. Reference Yu, Dong, Liu, Tang, Yuan and Xu2023). Recovery in the post-normal shock region of a single transonic interaction instead leads to outer layer growth according to the predicted adverse pressure gradient power law, as the flow adapts to the non-zero yet monotonically increasing outer pressure field (Pirozzoli et al. Reference Pirozzoli, Bernardini and Grasso2010). A more intricate path is displayed by the boundary layers in the multi-shock conical interaction of Zuo et al. (Reference Zuo, Memmolo, Huang and Pirozzoli2019), with multiple extended regions of APG and FPG acting on the outer boundary layer that lead to a gradual change between ZPG, APG and FPG states in the near-wall flow. The shock train case instead deviates from a near-ZPG state and follows an asymptotically decaying, oscillatory path in the phase space, while tending towards the theoretically predicted APG-type behaviour.

4.4. Momentum transport in the mean flow

Each interaction necessarily introduces changes in mean momentum transport as flux gradients develop in response to stresses imposed by outer shock waves and global boundary conditions. For compressible flows, the exact nature of these changes is assessed by decomposing the Favre-averaged streamwise momentum equation in the following manner (Anderson et al. Reference Anderson, Tannehill, Pletcher, Munipalli and Shankar2020):

(4.7) \begin{equation} \frac {\partial \bar {\rho }\tilde {u}}{\partial t} = - \bigg ( \frac {\partial \bar {\rho }\tilde {u}\tilde {u}}{\partial x} + \frac {\partial \bar {\rho }\tilde {u}\tilde {v}}{\partial y} \bigg ) -\frac {\partial \bar {p}}{\partial x} - \bigg ( \frac {\partial \overline {\rho }\widetilde {u^{\prime \prime }u^{\prime \prime }}}{\partial x} + \frac {\partial \overline {\rho }\widetilde {u^{\prime \prime }v^{\prime \prime }}}{\partial y} \bigg ) + \frac {1}{{\textit{Re}}_\infty } \bigg ( \frac {\partial \tilde {\tau }_{xx}}{\partial x} + \frac {\partial \tilde {\tau }{xy}}{\partial y} \bigg ), \end{equation}

allowing for budgets of $\rho \tilde {u}$ transport processes to be constructed. In (4.7), the terms on the right-hand side represent streamwise and vertical convective fluxes, inviscid pressure forces, apparent forces due to Reynolds stress variations and viscous diffusion effects. The left-hand side connotes the residual after summation over the right-hand side.

Distributions of these momentum budgets at the region bounds of table 1 are shown in figure 23. In the canonical boundary layer portion of the flow (figure 23 a), the prominent inner layer transport mechanisms are wall-normal viscous stresses and vertical gradients of turbulent shear stress, which are known to act against each other in wall-bounded turbulence (Jiménez Reference Jiménez2013). Secondary effects are present from convection, with wall-normal and streamwise gradients of momentum flux balancing one another as the boundary layer spatially develops along the channel wall. Across the inner and outer portions of the upstream boundary layer, effects of streamwise pressure gradient are negligible compared with the more dominant mean flow behaviours.

Substantial variations in momentum transport are generated in the flow after the leading shock boundary layer interaction in the channel. As shown in figure 23(b), strong streamwise pressure gradients are present throughout the local boundary layer. At the wall, these pressure forces are balanced by local vertical increases in molecular and turbulent shear. Beyond the near-wall region, turbulent shear continues to increase, enforcing net conservation of $\overline {\rho }\tilde {u}$ outside of the local viscous sublayer. As the wall-normal coordinate continues to increase, convective transport tends to dominate, which can be understood as the influx of high-momentum fluid introduced by a strong downward vertical velocity being balanced by local re-acceleration of the streamwise velocity (see figures 14, 15). Despite the strong shock-induced gradients, streamwise gradients in the Reynolds stress within the boundary layer are relatively small, possibly a consequence of longer time scale associated with adjustments in turbulence compared with the mean flow (Bevilaqua & Lykoudis Reference Bevilaqua and Lykoudis1978). Notable imbalances are observed in the outermost layer of the wall-bounded flow. However, this may be expected as the time-averaged LES, which is subjected to smearing from low frequency shock motions, is not strictly a converged fixed point solution to the RANS equations.

Figure 23. Streamwise momentum budgets for the shock train interaction at six representative streamwise stations. (a) $x/h=12$ . (b) $x/h=15.8$ . (c) $x/h=17.5$ . (d) $x/h=18.9$ . (e) $x/h=20$ . ( f) $x/h=25$ . Budgets are recorded from the wall to the local boundary layer edge. Negative values indicate terms that act as sinks on the mean $\overline {\rho }\tilde {u}$ , while positive values indicate source terms. As a result of boundary layer thickness fluctuations within the shock train, the presented budgets extend over various wall normal heights. To highlight differences between adjacent budgets, the vertical coordinate ( $y$ ) is normalised by the upstream viscous length $\ell _{v,{\boldsymbol{ref}}}$ .

Momentum budgets downstream of subsequent STBLIs in the shock train (figures 23 c–23 f) display similar distributions as those of figure 23(b), i.e. downstream STBLIs result in similar local mean flux gradients as the leading one. However, subtle differences emerge in the momentum transport mechanisms downstream as a result of the previously discussed gradual decay in pressure gradient imposed by the STBLIs as the local centreline Mach number approaches unity. Smaller streamline curvature is then sufficient to support wall-normal pressure gradients, which in turn weakens the outer layer momentum flux gradients due to convection. Similarly, outer layer variations in turbulent shear stress also diminish in response to the reduced outer pressure forces.

In the absence of strong outer layer perturbations, momentum transport is altered in the downstream region as the flow recovers in the weak adverse pressure gradient region. Decay of the balance term in figure 23( f), compared with upstream stations, indicates outer layer non-equilibrium in turbulent flow begins to dissipate and the strongest contributions to local momentum conservation once again reside within the inner layer. In this near-wall region, molecular and turbulent shear flux balance minor pressure gradient and convective effects. Increase in these viscous and buffer layer terms compared with prior stations is also attributed to the reintroduction of canonical near-wall dynamics, as discussed at length in the next section.

5.1. Reynolds stress evolution

In addition to altering the mean profiles of the primitive variables, the shock train causes significant modification to the turbulent stresses. Spatial distributions of all relevant Reynolds stress components are shown in figure 24. In each plot, the Reynolds stress components, given by the ensemble average of local density weighted velocity variances, are normalised by the square of the friction velocity at the upstream reference station, $x/h=12$ . All measured stresses undergo an increase in the vicinity of each STBLI, with the strongest increase occurring near the leading strong interaction. The chained line denotes the local maximum of turbulent kinetic energy, given by $({1}/{2})\widetilde {u_{i}^{\prime \prime }u_{i}^{\prime \prime }}$ , which is used as an approximate metric to identify the midpoint of the local shear layer.

Figure 24. Contour plots of the Reynolds stress components for the current spanwise periodic problem, normalised by the square of the reference station friction velocity. (a) Streamwise- normal Reynolds stress, $\widetilde {u^{\prime \prime }u^{\prime \prime }}/u_{\tau ,{\boldsymbol{ref}}}^2$ . (b) Wall- normal Reynolds stress, $\widetilde{\,v^{\prime\prime }v^{\prime\prime}}/u_{\tau ,{\boldsymbol{ref}}}^2$ (c) Spanwise- normal Reynolds stress, $\widetilde {w^{\prime \prime }w^{\prime \prime }}/u_{\tau ,{\boldsymbol{ref}}}^2$ . (d) Turbulent kinetic energy, $({1}/{2})\widetilde {u^{\prime \prime }_iu^{\prime \prime }_i}/u_{\tau ,{\boldsymbol{ref}}}^2$ . (e) Turbulent shear stress, $-\widetilde {u^{\prime \prime }v^{\prime \prime }}/u_{\tau ,{\boldsymbol{ref}}}^2$ .

The changes in Reynolds stress distributions are more quantitatively described by viewing the variation of each component in the wall-normal direction. Figure 25 presents such profiles at $10$ equally spaced streamwise locations over ${R}_1\to {R}_4$ . Each of the normalised $\widetilde{\,{u}^{\prime\prime}{u}^{\prime\prime}},\,\widetilde{\,{v}^{\prime \prime }{v}^{\prime \prime }},\,\widetilde {{w}^{\prime \prime }{w}^{\prime \prime }}$ and $\widetilde{\,{u}^{\prime\prime}{v}^{\prime\prime}}$ profiles retain expected TBL values and shapes upstream of the shock train. Strong amplification occurs as the leading oblique leg in the first STBLI is crossed; however, peaks of each stress component remain within the local inner layer ( $y/\delta \lt 0.1$ ). At downstream stations where the mean layer is thick and inflected, wall-normal locations of maximum stress shift to the local outer layer, a trend which holds for all stress tensor terms. Further downstream, the maxima in the profiles remain in the outer layer, though the amplitude of Reynolds stress terms begin to decay.

Figure 25. Profiles of the Reynolds stress tensor in the wall normal direction at collected streamwise stations. (a) $\widetilde {u^{\prime \prime }u^{\prime \prime }}^+$ . (b) $\widetilde{\,v^{\prime\prime }v^{\prime\prime}}^+$ . (c) $\widetilde {w^{\prime \prime }w^{\prime \prime }}^+$ . (d) $-\widetilde {u^{\prime \prime }v^{\prime \prime }}^+$ . Vertical locations are normalised using the local $\delta$ , while stress values are normalised by $u_{\tau ,{\boldsymbol{ref}}}^2$ . This approach highlights both the shift of stresses from the inner layer to the outer layer and the absolute changes in stress magnitude.

5.2. Turbulence amplification

A quantitative description of turbulence amplification further clarifies the influence of the shock train on the incoming boundary layer. For a given generic stress component $\widetilde {u_{i}^{\prime \prime }u_{j}^{\prime \prime }}$ , an amplification factor is defined as

(5.1) \begin{equation} \varLambda _{\textit{ij}}(x) = \frac {|\widetilde {u_{i}^{\prime \prime }u_{j}^{\prime \prime }}(x,y)|_\infty }{|\widetilde {u_{i}^{\prime \prime }u_{j}^{\prime \prime }}(x_0,y)|_\infty }, \end{equation}

where $x_0$ is a reference streamwise location and $|\boldsymbol{\cdot }|_\infty$ is the usual max norm.

Figure 26(a) displays the amplification of all individual Reynolds stress components and the TKE, taking the reference location $x_0$ to be in the undisturbed boundary layer at the end of ${R}_1\,(x/h=12)$ . Strong amplification is observed beneath the leading STBLI for all components, with the wall-normal and spanwise components being augmented the most. Peak amplification occurs for all components in this region except $\widetilde{\,{u}^{\prime\prime}{v}^{\prime\prime}}$ , which attains its local maxima within the interaction between the second shock and the thickened wall shear layer, as also noted by Carroll & Dutton (Reference Carroll and Dutton1992). Furthermore, each component exhibits an oscillatory amplification and attenuation pattern as individual STBLIs are traversed. The $\widetilde{\,{u}^{\prime\prime}{u}^{\prime\prime}},\,\widetilde {{w}^{\prime \prime }{w}^{\prime \prime }},\,\widetilde{\,{u}^{\prime\prime}{v}^{\prime\prime}}$ and TKE curves attain local maxima at approximately the midpoint of each interaction region. The $\widetilde{\,{v}^{\prime \prime }{v}^{\prime \prime }}$ curves display different behaviour, instead attaining local minima near each midpoint. Both the streamwise stress and the TKE, after being initially amplified, decay below their reference value at far downstream stations. The spanwise, wall-normal and shearing components remain amplified above their reference values. A summary of the amplification behaviour for each stress tensor component is provided in table 4.

Figure 26. Turbulence amplification characteristics shown by (a) the maximum amplification value at each streamwise station and (b) the wall-normal location of each streamwise maximum.

Table 4. Maximum amplification factors and locations for the normal Reynolds stresses, the turbulent kinetic energy and the turbulent shear stress.

The wall-normal locations $(y/h)$ of each value of $\varLambda _{\textit{ij}}(x)$ are shown in figure 26(b). In the upstream region, TKE, turbulent shear and wall parallel stress components are strongest near the wall. The wall-normal component is largest slightly farther from the wall, owing to kinematic constraints introduced by the no penetration condition at the boundary. After the first STBLI is crossed, all components are strongest far away from the wall, as the shear layer thickens under the influence of the shock train. Returning to the distributions in figure 25, these locations are readily identified in the outer region of the downstream, recovering boundary layer. The physical height of maximum amplification of the Reynolds shear stress deviates from the other components in the vicinity of reattachment, at $x/h\sim 15.5$ , indicating a slightly varied response in the turbulent shear stress compared with the normal stresses.

5.3. Stress tensor anisotropy

One measure of turbulence alterations through the sequence of STBLI is based on the underlying structure of the Reynolds stress tensor, whose components are differentially amplified. Rigorous investigation of these changes is undertaken by using Lumley triangles, offering a point-wise examination of anisotropy that complements the purely statistical analysis from the previous sections. Anisotropy maps are computed at six representative streamwise stations, given by the spatial bounds (see table 1) of ${R}_1,\,{R}_2,\,{R}_3,\,{R}_4\,$ and the downstream region. Normalised anisotropic stress components are derived from the original Reynolds stress tensor (Pope Reference Pope2001) according to

(5.2) \begin{gather} b_{\textit{ij}} = \frac {R_{\textit{ij}}}{\textit{trace}(R_{\textit{ij}})} - \frac {1}{3}\delta _{\textit{ij}} = \frac {\widetilde {u^{\prime \prime }_iu^{\prime \prime }_j}}{2k} - \frac {1}{3}\delta _{\textit{ij}}. \end{gather}

By computing the two non-zero invariants ( $I_2,\,I_3$ ) associated with the second-order tensor $b_{\textit{ij}}$ , intrinsic directional preference to the turbulent stresses may be educed. All realisable combinations of $I_2$ and $I_3$ fall within an approximately triangular shape. Under the adopted invariant definitions, the turbulent fluctuations have the limiting states of $(I_2,I_3) = (0,0)$ , $(I_2,I_3) = (({-1}/{3}),({2}/{27}))$ and $(I_2,I_3) = (({-1}/{12}),({-1}/{108}))$ . These respectively correspond to isotropic, one component and two component axisymmetric states of the Reynolds stress tensor (Simonsen & Krogstad Reference Simonsen and Krogstad2005). Note that when plotting, $-I_2$ is used in place of $I_2$ for visualisation convenience.

Wall-normal distributions of these invariants, shown in figure 27, quantitatively reveal a reorganisation of the near-wall turbulent flow across each region of the shock train. Upstream of the interaction in ${R}_1$ (figure 27 a), the flow displays a purely two-component preference immediately at the wall, before tending towards a one-component state within the local buffer region. Such behaviour has been observed in canonical wall-bounded flows, where in the viscous sublayer, wall-normal fluctuations are greatly damped and in the buffer layer, the autonomous streak cycle generates dominant streamwise aligned velocity disturbances (Jiménez & Moin Reference Jiménez and Moin1991). Moving away from the wall, the turbulence field becomes increasingly isotropic, as both the kinematic damping of the wall and the mean shear decrease.

Downstream of the first STBLI at the upper bound of ${R}_1$ (figure 27 b), the turbulence exhibits a two component axisymmetric state near the wall. The structure of the invariant map indicates that while wall-normal fluctuations remain suppressed, the near-wall streak cycle has been altered such that streamwise and spanwise fluctuations are now of similar magnitude. Confirmation of this observation is seen in figure 25, as downstream of the first STBLI, $\widetilde{\,{u}^{\prime\prime}{u}^{\prime\prime}}$ and $\widetilde {{w}^{\prime \prime }{w}^{\prime \prime }}$ fluctuations display similar amplitudes, while the $\widetilde{\,{v}^{\prime \prime }{v}^{\prime \prime }}$ stresses remain comparably weak in the local inner layer $(y\leqslant 0.1\delta )$ . The first STBLI also alters the structure of the outer layer, with a more isotropic state being measured at the end of ${R}_1$ than in the undisturbed upstream turbulence.

Invariant maps downstream of the second (figure 27 c), third (figure 27 d) and fourth (figure 27 e) STBLIs show slightly different inner layer structures than those downstream of the leading STBLI. While a clear two-component preference is maintained close to the wall, a modest re-emergence of a one-component preference is evident in the increasingly bulging profile of the invariant map in the inner layer. Outer layer stress isotropy is maintained as the thick wall layer interacts with these downstream STBLIs, indicating the adverse pressure gradient imposed by the shock train generates and sustains an equipartition of the turbulent stresses far from the wall.

At the most downstream station (figure 27 f), the tendency towards single-component anisotropy becomes exaggerated in the inner layer. The increase in anisotropy with streamwise distance supports the hypothesis that as the shock train is traversed, the flow patterns of the inner layer recover gradually towards those of a canonical wall-bounded flow. Reorganisation in the downstream region is further facilitated by the weakening of the local adverse pressure gradient (figure 5), as the STBLIs in this portion of the flow occur at near sonic Mach numbers (figure 12).

Figure 27. Normalised anisotropy maps of the Reynolds stress tensor, at six representative streamwise stations. (a) $x/h=12$ . (b) $x/h=15.8$ . (c) $x/h=17.5$ . (d) $x/h=18.9$ . (e) $x/h=20$ . ( f) $x/h=25$ . Insets are included to give a clearer picture of the return to isotropy near the origin of the invariant maps.

Redistribution of the stress components from an anisotropic state to an isotropic one is a statistical result of a structural change in wall turbulence dynamics. Figure 28 demonstrates this process, presenting the streamwise velocity field in an $x{-}z$ plane at a height of $y^+_{\textit{ref}} = 12$ . Low-speed streaks, characteristic of an autonomous wall cycle (Jiménez & Moin Reference Jiménez and Moin1991), are identified in the upstream region. These structures leave the plane as the flow separates due to the leading STBLI, with near-wall flow in the vicinity of the weaker STBLIs also intermittently interrupted by flow separation. However, as the interaction strength decreases (figure 18), the near-wall flow is increasingly unperturbed by inviscid shock waves and the mean shear begins to recover (figure 17). At far-downstream stations, where the flow is largely unaffected by shock-induced pressure gradients, streak-like structures begin to re-emerge in the flow.

Figure 28. Near-wall view of instantaneous streamwise velocity field at a wall-normal height corresponding to $y^+_{\textit{ref}}\sim 12$ . Shown are values of $u/u_\infty$ ranging from $0\, \mathrm{to}\, 0.6$ . Regions of instantaneous reverse flow are outlined in red, and white vertical lines denote individual STBLI region boundaries.

5.4. Turbulent kinetic energy budgets

Further assessment of the effect of shock trains on the turbulent near-wall flow is obtained by examining TKE budgets at various points. Budgets highlight the complex transport processes that emerge as the turbulence is repeatedly altered within the shock train, which can help furnish turbulence models for lower fidelity approaches such as RANS or simplified analytical solutions (Babinsky & Harvey Reference Babinsky and Harvey2011). The transport equation for time-averaged turbulent kinetic energy is

(5.3) \begin{equation} \frac {\partial \bar {\rho }\tilde {k}}{\partial t} = C + T + P + V -\bar {\rho }\epsilon + K, \end{equation}

with $C$ , $T$ , $P$ , $V$ , $-\bar {\rho }\epsilon$ and $K$ representing mean flow convection, turbulent transport, production, viscous diffusion, turbulent dissipation and genuine compressibility effects. The form of each individual term is given according to Pirozzoli et al. (Reference Pirozzoli, Bernardini and Grasso2010) as

(5.4) \begin{align} C &= -\frac {\partial \bar {\rho }\tilde {u}_j \tilde {k}}{\partial x_j }, \end{align}

(5.5) \begin{align} T &= -\frac {\partial }{\partial x_j } \bigg ( \frac {1}{2}\widetilde {u^{\prime \prime }_iu^{\prime \prime }_iu^{\prime \prime }_j} + \overline {p'u^{\prime \prime }_j} \bigg ), \end{align}

(5.6) \begin{align} P &= -\overline {\rho }\widetilde {u_{i}^{\prime \prime }u_{j}^{\prime \prime }}\frac {\partial \tilde {u}_i}{\partial x_j}, \end{align}

(5.7) \begin{align} V &= \frac {\partial }{\partial x_j} (\overline {\sigma _{\textit{ij}}'u^{\prime \prime }_i}), \end{align}

(5.8) \begin{align} \bar {\rho }\epsilon &= \overline {\sigma _{\textit{ij}}' \frac {\partial u^{\prime \prime }_i}{\partial x_j} } ,\end{align}

(5.9) \begin{align} K &= \overline {p'\frac {\partial u^{\prime \prime }_i}{\partial x_i}} + \overline {u^{\prime \prime }_i} \bigg ( \frac {\partial \overline {\sigma _{\textit{ij}}}}{\partial x_j} - \frac {\partial \overline {p}}{\partial x_j} \bigg ). \end{align}

Figure 29. Turbulent kinetic energy budgets for the shock train interaction at six representative streamwise stations. (a) $x/h=12$ . (b) $x/h=15.8$ . (c) $x/h=17.5$ . (d) $x/h=18.9$ . (e) $x/h=20$ . ( f) $x/h=25$ . Note the varying abscissa scales, with panel (a) displaying only near-wall budget variations, as this is the streamwise station where equilibrium TBL behaviour is observed. Other stations display the entirety of the budget across the recorded local boundary layer height. Profiles in panel (a) are compared with DNS data (coloured squares) from Pirozzoli et al. (Reference Pirozzoli, Grasso and Gatski2004).

Normalised budgets are shown in figure 29. At the most upstream station (figure 29 a), only near-wall data are presented as most TKE transport occurs within the inner layer (Pope Reference Pope2001). Favourable agreement is observed between phenomena of the current low-Re TBL and the equilibrium Mach 2.25 TBL of Pirozzoli, Grasso & Gatski (Reference Pirozzoli, Grasso and Gatski2004). Viscous diffusion and dissipation terms attain their global extrema at the wall, and decay significantly throughout the linear sublayer ( $y^+\leqslant 5$ ). In the buffer layer, the production term attains its maximum value, and the local gain in TKE is balanced by viscous diffusion, turbulent transport and turbulent dissipation. In the log and outer layers ( $y^+\gtrapprox 30$ ), only production and dissipation are present, with equilibrium being ensured by $P/\overline {\rho }\epsilon \sim 1$ . The effects of genuine compressibility and mean flow convection are negligible throughout the undisturbed boundary layer.

After the leading bifurcated STBLI, the distributions of the budget terms are altered significantly. Most notably, TKE production now occurs in the outer layer of the flow. Turbulent dissipation nonetheless balances production in this region, but in contrast to the upstream TBL, production in the outer layer is also counteracted by mean convection and turbulent transport. Budget terms in the inner layer demonstrate the influx of TKE through both viscous diffusion and turbulent transport is locally balanced by dissipation, with contributions from local production being largely muted in this region.

Budgets at the downstream ends of ${R}_2$ , ${R}_3$ , ${R}_4$ and beyond bear similarity to the budget recorded after the leading STBLI. However, multiple trends develop in both the inner and outer layers over the streamwise distance spanned by the budget locations. Production in the outer layer decreases significantly along the streamwise direction. To maintain outer layer balance, the effects of mean convection, turbulent transport and dissipation also weaken. Within the inner layer, a secondary peak in production emerges as the streamwise coordinate increases, with local balance enforced by molecular diffusion, dissipation and turbulent fluctuations. This double hump structure to the production profile may be viewed as evidence for recovery of the wall streak cycle, corroborated by the preceding anisotropy maps and near-wall visualisation.

5.5. Production of turbulent kinetic energy

Turbulent kinetic energy budgets and anisotropy maps show that the processes by which turbulent fluctuations are produced varies as a function of both $x$ and $y$ within the shock train. To clarify this variation, the production term from the TKE transport equation,

(5.10) \begin{equation} {P}_k = -\overline {\rho }\widetilde {u_{i}^{\prime \prime }u_{j}^{\prime \prime }}\frac {\partial \tilde {u}_i}{\partial x_j}, \end{equation}

is further analysed. For spanwise homogeneous flows, total production can be further decomposed into its streamwise, wall-normal and shearing components according to

(5.11) \begin{equation} {P}_k = {P}_x + {P}_y + {P}_s = -\overline {\rho }\widetilde{\,{u}^{\prime\prime}{u}^{\prime\prime}}\frac {\partial \tilde {u}}{\partial x} -\overline {\rho }\widetilde {{v}^{\prime \prime }{v}^{\prime \prime }\ }\frac {\partial \tilde {v}}{\partial y} -\overline {\rho }\widetilde {{u}^{\prime \prime }{v}^{\prime \prime }\ }\bigg (\frac {\partial \tilde {u}}{\partial y} + \frac {\partial \tilde {v}}{\partial x}\bigg ). \end{equation}

The above-mentioned partitioning is invoked to evaluate the maximum values and wall-normal locations of ${P}_k, {P}_x, {P}_y$ and ${P}_s$ at each streamwise location. Analysis of these production terms is restricted to points within the local boundary layer, i.e. ${P}(x,y) \in \{ (x,y) : y\leqslant \delta (x) \}$ , to avoid computing spurious ‘peaks’ in production that are a statistical artefact of shock wave oscillations rather than true turbulence amplification mechanisms.

Figure 30. Total and individual production characteristics for TKE shown by (a) the maximum value of each production mechanism attained within the boundary layer at each streamwise station and (b) the wall-normal location of each streamwise maximum.

Figure 30(a) demonstrates that in the upstream boundary layer, TKE is produced entirely through shearing mechanisms, associated with the well-described sweep/ejection events that dominate undisturbed wall bounded shear flows (Jiménez Reference Jiménez2013). As the leading viscous separation region ( $x/h = 13.72$ ) is approached, total production increases rapidly, marked by an initial peak in ${P}_s$ followed immediately downstream by the apex of the ${P}_x$ profile. Similar results were observed in the impinging oblique STBLI of Fang et al. (Reference Fang, Zheltovodov, Yao, Moulinec and Emerson2020), where production due to streamwise stress-mean flow interactions was the dominant TKE source over the leading portion of the separated region. Owing to wall-normal deflection of the boundary layer and area confinement over the current separation bubble, the vertical production term acts as a non-negligible sink near the separated flow region.

Near the second STBLI, total TKE production displays a local maximum as the turbulent stresses act against the shock induced pressure gradients. In contrast to the first STBLI, this peak in total production comprises largely a shear-type mechanism, indicating subtle differences in the turbulent response between the foremost and secondary STBLIs in the train. However, a spatial distribution of production variables similar to the leading STBLI is observed, with a crest in ${P}_s$ slightly preceding a crest in ${P}_k$ , which is rapidly followed by a synchronised crest/trough pair in ${P}_x$ and ${P}_y$ underneath the outer layer shock wave. This pattern repeats for all subsequent STBLI, with production amplitude variations decreasing in magnitude as the interactions become weaker. Beyond $x/h\sim 23$ , production is comprised almost entirely by the shear term, ${P}_s$ , indicating production behaviours in the recovering boundary layer share commonalities with those in the undisturbed boundary layer upstream of ${R}_1$ .

The wall-normal location of each production term maximum varies considerably over the shock train interaction. As seen in figure 30(b), after deflection of the boundary layer in ${R}_1$ , all production terms retain peaks in outer layer, though each mechanism displays a slightly different $y(x)$ variation. Shifts of large production values to the local outer layer are in agreement with the finding that downstream of separation, strong turbulent fluctuations reside in the outer portions of the thickened wall boundary layer. Another notable trend from figure 30(b) is that shear and total production mechanisms are nearly coincident in the outer layer from the point of separation until $x/h\sim 24$ , whereupon a dramatic decrease in the height of peak production is observed. This large wall-normal shift in the location of peak TKE production arises from the combination of gradual growth of near-wall shear and a gradual decay of outer layer turbulent stresses, previously noted in figure 26. Maximisation of the shear and total production terms near the wall, given the mean shear is positive in this region, indicates a strong, negative $\widetilde {{u}^{\prime \prime }{v}^{\prime \prime }\ }$ correlation is developing close to the wall that is responsible for the local production of TKE.

Dominance of shear-induced production, considered in conjunction with the already discussed stress anisotropy, TKE budget and skin friction recovery results, provide evidence for the re-emergence of the classical near-wall streak cycle. However, the Reynolds stress amplification analysis indicates that in the region of $x/h\gt 24$ , the strongest turbulence stresses still exist well within the local outer layer. The disparity between production maximising near the wall and Reynolds stresses maximising in the outer layer is likely a turbulence memory effect, where the outer regions of the flow inherit larger energetic scales from upstream stations while the re-energised, small near-wall scales in the weak local pressure gradient are expected to re-establish canonical inner and overlap layer behaviours as they propagate downstream.

Figure 31. Pre-multiplied spanwise energy wavelength spectra as a function of wall-normal distance. (a) $x/h=12$ . (b) $x/h=15.8$ . (c) $x/h=17.5$ . (d) $x/h=18.9$ . (e) $x/h=20$ . ( f) $x/h=25$ . The abscissa uses the local viscous length scale from the upstream boundary layer ( $x/h=12$ ). Spectra are not normalised and, as such, are presented in arbitrary units based on local spectral magnitudes.

Figure 32. Pre-multiplied streamwise energy wavelength spectra as a function of wall-normal distance. (a) $x/h=12$ . (b) $x/h=15.8$ . (c) $x/h=17.5$ . (d) $x/h=18.9$ . (e) $x/h=20$ . ( f) $x/h=25$ . The abscissa uses the local viscous length scale from the upstream boundary layer ( $x/h=12$ ). Spectra are not normalised and, as such, are presented in arbitrary units based on local spectral magnitudes.

5.6. Modification of length scales

Statistical changes in the wall turbulence and the emergence of varied inner and outer layer processes are associated with physical changes to the boundary layer structures. This notion is explored by computing the pre-multiplied spanwise ( $\varPhi _{uu} = k_z E_{uu}(k_z)$ ) and streamwise ( $\varPhi _{uu} = k_x E_{uu}(k_x)$ ) spectral densities (Smits et al. Reference Smits, McKeon and Marusic2011) of the streamwise velocity variances, allowing a quantification of the most energetic length scales. Results from six streamwise locations are compared, corresponding to the bounds of ${R}_1\to {R}_4$ and an additional station at the end of the downstream region. The resulting distributions are plotted in figures 31 and 32, where the spanwise spectra are computed by invoking spanwise periodicity and the streamwise spectra computed using the Taylor hypothesis. Locations of the most energetic scales, and their wall-normal locations, are given in table 5. Units are normalised by the viscous length scale from the upstream boundary layer, rendering the spectra functions of scaled wavelength ( $\lambda _{x_i}^+ = 1/k_{x_i}^+$ ) and scaled height ( $y^+_{\textit{ref}}$ ).

In the undisturbed TBL upstream of the shock train, the most energetic spanwise wavelengths are $\lambda _z^+\sim 100$ viscous lengths and the most energetic streamwise wavelengths are $\lambda _x^+\sim 1500$ viscous lengths. Both spectral peaks occur within the local buffer layer, at $y^+\sim 10$ . These physics are well aligned with those documented for equilibrium turbulent boundary layers, where the process of near-wall streak formation and vortex lifting contributes most strongly to the streamwise velocity variance (Perry & Chong Reference Perry and Chong1982). Note in figures 31(a) and 32(a), only inner peaks are observed as a consequence of the considered low Reynolds number.

Downstream of the leading STBLI, an appreciable shift appears in the location of both the spanwise and streamwise energetic wavelengths. The most energetic scales now reside in the outer layer of the flow, at a height of $y^+_{\textit{ref}}\sim 200$ . Despite the location shift, energy remains concentrated at the same streamwise scales measured in the upstream TBL region, with peak variance residing near $\lambda ^+_x\sim 1000$ . Energy in the spanwise scales however resides in wavelengths near $\lambda ^+_z\sim 250$ , indicating that while the streamwise distribution of structures is inherited over the leading interaction zone, a spanwise pairing of adjacent eddies may be present in this region.

Table 5. Size and wall-normal location of most energetic length scales throughout the shock train interaction. All scales have been normalised using wall units from the upstream TBL station.

Despite interactions with subsequent shock waves, no abrupt changes in the energy wavelength spectra are observed at the ends of ${R}_2$ , ${R}_3$ , ${R}_4$ , or the downstream region. Instead, a gradual increase in both the streamwise and spanwise energetic length scales is observed. The spanwise wavelength spectrum at the end of the downstream region, shown in figure 31( f), has an energetic peak at $\lambda _z^+\sim 480$ , based on the upstream reference viscous scales. Similarly, the streamwise wavelength spectrum of figure 32( f) has its peak at $\lambda _x^+\sim 3900$ , indicating energy is shifted into scales approximately $3$ to $4$ times larger than those in the undisturbed boundary layer.

A continual increase in the wall-normal location of peak variance is also observed in the energy wavelength spectra. In correspondence with the increase in the local boundary layer thickness, the vertical location of these energetic scales shifts slightly outward, moving from $y^+_{\textit{ref}}\sim 200$ after the leading STBLI to $y^+_{\textit{ref}}\sim 385$ . The wall-normal movement of this spectral peak indicates that the most energetic turbulent fluctuations occur within the outer layer, consistent with the results from previous sections focusing on Reynolds stress amplification and TKE budgets. Further evidence for the recovery of the near-wall turbulence can be seen in figure 31( f), where a secondary peak at $(y^+_{\textit{ref}},\lambda ^+_z) \sim (10,100)$ re-emerges that corresponds to the repeated formation of low-speed streaks in the local buffer layer (Yu et al. Reference Yu, Dong, Liu, Tang, Yuan and Xu2023).

There are also non-local, or imprinting, artefacts observed by the large outer layer eddies on the resurgent near-wall cycle. In figures 31( f) and 32( f), spectral handles can be seen to reach from the energetic larger scales to the local buffer layer. As the buffer layer cycle exists autonomously, individual bursting events in this region can effectively reside within the footprints of large scales, provided a certain scale separation exists. The natural generation of larger structures by the shock train interaction provides such a mechanism, and as these eddies persist downstream (figure 32) smaller near-wall features may persist within them.

We conclude the analysis of the near-wall turbulence in a shock train by presenting a qualitative visualisation of the fluctuating streamwise velocity field. Slices from an instantaneous realisation of the simulation are shown at increasing wall-normal distances in figure 33, and this figure serves to convey in a more intuitive manner the physics evidenced throughout § 5. In figure 33(a), smaller scale streaks are clearly present upstream of the shock train interaction. These are observed to lift away from the wall over the leading separation region (figure 33 b), with certain aspects of the fine scale organisation being inherited by downstream locations (figure 33 c). Beyond the lift up location, however, the fine grain coherence of the turbulent structures begins to rapidly decay and is replaced by qualitatively larger structures. These larger scales persist for a majority of the downstream portion of the shock train and exist well within the local outer layer (figures 33 d– 33 e). However, at far downstream regions, near-wall streaks can even visually be seen to reappear in the inner layer, commensurate with the downstream region where flow recovery initiates (figure 33 a).

Figure 33. Visualisations of instantaneous streamwise velocity perturbations ( $u'/u_\infty$ ) at ascending wall-normal heights. (a) $y/h = 0.01, y^+_{\textit{ref}} = 12$ . (b) $y/h = 0.09, y^+_{\textit{ref}} = 80$ . (c) $y/h = 0.17, y^+_{\textit{ref}} = 158$ . (d) $y/h = 0.26, y^+_{\textit{ref}} = 237$ . (e) $y/h = 0.41, y^+_{\textit{ref}} = 372$ . Vertical red lines demarcate region boundaries associated with individual STBLIs in the flow.