3.1. Mathematical formulation

3.1.1. Preshock perturbations as a composition of mono-frequency modes

The formulation employed in the present study is largely consistent with that of Wouchuk et al. (Reference Wouchuk, Huete Ruiz de Lira and Velikovich2009), Huete et al. (Reference Huete, Cuadra, Vera and Urzay2021), Cuadra et al. (Reference Cuadra, Vera, Di Renzo and Huete2023) and Cuadra (Reference Cuadra2023), in which the upstream flow consisted solely of vortical disturbances. Here, the model is extended to account for density (Huete et al. Reference Huete, Velikovich and Wouchuk2011) and pressure (Huete et al. Reference Huete, Wouchuk and Velikovich2012b ) fluctuations in the upstream flow, associated with acoustic, entropic and vortical disturbances. Furthermore, the effect of vibrational excitation is incorporated through the solution of the base-flow equations in (2.1), using the Combustion Toolbox (Cuadra et al. Reference Cuadra, Huete and Vera2025b , Reference Cuadra, Huete and Vera2026), a general-purpose framework for thermochemical modelling of multi-component mixtures. Among its capabilities, the code can solve problems involving reacting and non-reacting shocks under various gas models. In this work, we adopt a thermally perfect gas model in which vibrational excitation is captured through the NASA-9 polynomial fits (McBride Reference McBride2002) embedded in the toolbox. For a broader analysis in terms of shock intensity, we refer to Cuadra et al. (Reference Cuadra, Di Renzo, Hoste, Williams, Vera and Huete2025a ).

The correlation between the nature of entropic and vorticity disturbances suggests the decomposition of the upstream turbulent flow into two types of waves: the vorticity–entropy wave and the acoustic wave. The former can be expressed as a steady wave function in the absence of diffusive effects, namely

(3.1a) \begin{align} \delta u_1^r \left ( x, y \right )&= \varepsilon _r \langle c_1\rangle \cos \big (k_x^r x \big )\cos \big (k_y y\big ), \\[-12pt]\nonumber \end{align}

(3.1b) \begin{align} \delta v_1^r \left ( x, y \right )&= \varepsilon _r \langle c_1\rangle \dfrac {k_x^r}{k_y} \sin \big (k_x^r x \big )\sin \big (k_y y\big ), \\[-12pt]\nonumber \end{align}

(3.1c) \begin{align} \delta \rho _1^e \left ( x, y \right ) &= \varepsilon _e \langle \rho _1\rangle \cos \big (k_x^r x \big )\cos \big (k_y y\big ); \end{align}

upon application of the divergence-free condition in the velocity field. The acoustic mode, corresponding to isentropic curl-free travelling waves, involves perturbations in the form

(3.2a) \begin{align} \delta u_1^a \left ( x, y,t \right )&= \varepsilon _a \langle c_1\rangle \cos \left (k_x^a x - \omega _a t\right )\cos \big (k_y y\big ) , \\[-12pt]\nonumber \end{align}

(3.2b) \begin{align} \delta v_1^a \left ( x, y,t \right )&= -\varepsilon _a \langle c_1\rangle \dfrac {k_y}{k_x^a} \sin \left (k_x^a x - \omega _a t\right )\sin \big (k_y y\big ), \\[-12pt]\nonumber \end{align}

(3.2c) \begin{align} \delta \rho _1^a \left ( x, y ,t\right ) &= \varepsilon _a \langle \rho _1\rangle \dfrac {\omega _a}{k_x^a} \cos \left (k_x^a x-\omega _a t \right )\cos \big (k_y y\big ) , \\[-12pt]\nonumber \end{align}

(3.2d) \begin{align} \delta p_1 \left ( x, y,t \right ) &= \varepsilon _a \langle \rho _1 \rangle \langle c_1 \rangle ^2 \dfrac {\omega _a}{k_x^a} \cos \left (k_x^a x-\omega _a t \right )\cos \big (k_y y\big ) . \end{align}

In this formulation, $\langle c_1\rangle$ and $\langle \rho _1\rangle$ denote the mean speed of sound and density in the preshock gas, corresponding to the perturbation-free conditions. The dimensionless parameter $\varepsilon _r$ denotes the preshock streamwise velocity rotational fluctuation amplitude, assumed to be small in the linear theory, $\varepsilon _r\ll 1$ . Characterised by a smaller level of intensity, the density-entropic fluctuation amplitude $\varepsilon _e$ and the streamwise velocity-acoustic fluctuation amplitude $\varepsilon _a$ are so defined. In Huete et al. (Reference Huete, Wouchuk, Canaud and Velikovich2012a ), the case of the lagged reshock problem – a second time-retarded shock interacting with the turbulence generated by a first shock moving in an entropic turbulent flow – the relation between the rotational and entropic disturbances was determined by the first shock properties. In the case considered in this work, the upstream turbulent flow is dominated by the rotational disturbances, so that $\varepsilon _e$ and $\varepsilon _a$ are smaller or much smaller than $\varepsilon _r$ . However, their impact on the postshock turbulent flow is not necessarily negligible, as will be demonstrated later. In real applications, such interrelation between the modes comes from additional effects that have not been considered in the LIA description, such as weak nonlinearities and dissipation effects. Therefore, the relation between the three types of perturbation amplitudes must be externally provided. The numerical turbulence generator, employed later in the DNS study, will be used to compute the relative amplitude of the respective modes ahead of the shock.

The rotational-entropic perturbation is also characterised by the longitudinal wavenumber $k_x^r$ , which is employed to define the upstream vorticity field, namely

(3.3) \begin{equation} \varpi _{z,1} \left ( x, y \right ) =\varepsilon _r \langle c_1\rangle k_y\left [1+\left (\frac {k_x^r}{k_y}\right )^2 \right ] \cos \big (k_x^r x \big )\sin \big (k_y y\big ). \end{equation}

The acoustic field demands the definition of the corresponding wavenumber $k_x^a$ and frequency $\omega _a$ , which are related to each other through the sonic dispersion relationship $\omega _a^2 = c_1^2[ (k_x^a)^2+k_y^2]$ . It is found convenient to define the wave angles $k_x^r/k_y=1/|\!\tan \theta ^r|$ and $k_x^a/k_y=1/|\!\tan \theta ^a|$ , which are used later in the analytical computation of the postshock mean values through the corresponding probability density functions. Note that a similar transverse wavenumber $k_y$ is used in both rotational-entropic and acoustic modes. This anticipated condition arises because all types of disturbances are generated through the same source and, therefore, share the main properties of the turbulent spectrum.

To perform the normal mode analysis, it is also found useful to write the upstream perturbation variables in the shock reference frame. For this purpose, let us first define the dimensionless time as $\tau = t k_y c_2$ , so the upstream excitation ahead of the shock translates into constant-amplitude oscillations in the asymptotic long-time regime, where

(3.4) \begin{equation} \omega _{s}^r= \dfrac {{\mathcal R} {\mathcal M}_2}{\tan \theta ^r} \quad {\text{and}}\quad \omega _{s}^a= \dfrac {{\mathcal R} {\mathcal M}_2}{\tan \theta ^a}\left (1-\dfrac {1}{{\mathcal M}_1\cos \theta ^a}\right ) \end{equation}

refer to the dimensionless oscillation frequencies produced by the upstream rotational-entropic and acoustic modes, respectively. Thus, velocity, density, pressure and vorticity perturbations right ahead of the shock, in terms of the dimensionless time $\tau$ , take the form

(3.5a) \begin{align} \overline {u}_{1s}(\tau )&=\dfrac {\delta u_{1s}}{\varepsilon _r c_1}= \cos {\big(\omega _s^r\tau +\phi _r\big)}+ \dfrac {\varepsilon _a}{\varepsilon _r} \cos {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.5b) \begin{align} \overline {v}_{1s}(\tau )&=\dfrac {\delta v_{1s}}{\varepsilon _r c_1}= \dfrac {k_x^r}{k_y}\sin {\big(\omega _s^r\tau +\phi _r\big)}- \dfrac {k_y}{k_x^a}\dfrac {\varepsilon _a}{\varepsilon _r} \sin {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.5c) \begin{align} \overline {\rho }_{1s}(\tau )&=\dfrac {\delta \rho _{1s}}{\varepsilon _r \rho _1}= \dfrac {\varepsilon _e}{\varepsilon _r}\cos {\big(\omega _s^r\tau +\phi _e\big)}+ \dfrac {\omega _a}{k_x^a}\dfrac {\varepsilon _a}{\varepsilon _r} \cos {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.5d) \begin{align} \overline {p}_{1s}(\tau )&=\dfrac {\delta p_{1s}}{\varepsilon _r \rho _1 c_1^2}= \dfrac {\omega _a}{k_x^a}\dfrac {\varepsilon _a}{\varepsilon _r} \cos {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.5e) \begin{align} \overline {\varpi }_{1s}(\tau )&=\dfrac {\varpi _{1s}}{\varepsilon _r k_y c_1}= \left [1+\left (\frac {k_x^r}{k_y} \right )^2\right ]\cos {\big(\omega _s^r\tau +\phi _r\big)}, \end{align}

where the average symbol brackets $\langle \boldsymbol{\cdot }\rangle$ have been omitted for clarity. In writing (3.5), $\phi _r$ , $\phi _e$ and $\phi _a$ denote the phases of the rotational, entropic and acoustic modes, respectively, and the subscript $1s$ defines upstream (subscript $1$ ) disturbances at the shock (subscript $s$ ) location.

3.1.2. Postshock perturbations for dual-frequency composition modes

Once the upstream conditions are determined, the next step is to obtain the postshock properties for the different perturbation modes. In a system of reference moving with the shock front immediately downstream from the perturbed shock surface (subscript $2$ ), a corresponding set of perturbations is generated

(3.6a) \begin{align} \overline {u}_{2s}(\tau )&=\dfrac {\delta u_{2s}}{\varepsilon _r c_2}= \big({\mathcal U}_r^r+{\mathcal U}_a^r\big) \cos {\big(\omega _s^r\tau +\phi _r\big)}+ \big({\mathcal U}_r^a+{\mathcal U}_a^a\big)\cos {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.6b) \begin{align} \overline {v}_{2s}(\tau )&=\dfrac {\delta v_{2s}}{\varepsilon _r c_2}=\big({\mathcal V}_r^r+{\mathcal V}_a^r\big)\sin {\big(\omega _s^r\tau +\phi _r\big)} + \big({\mathcal V}_r^a+{\mathcal V}_a^a\big) \sin {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.6c) \begin{align} \overline {\rho }_{2s}(\tau )&=\dfrac {\delta \rho _{2s}}{\varepsilon _r \rho _2}=\big({\mathcal D}_e^r+{\mathcal D}_a^r\big)\cos {\big(\omega _s^r\tau +\phi _e\big)}+ \big({\mathcal D}_e^a+{\mathcal D}_a^a\big)\cos {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.6d) \begin{align} \overline {p}_{2s}(\tau )&=\dfrac {\delta p_{2s}}{\varepsilon _r \rho _2 c_2^2}={\mathcal P}_a^r\cos {\big(\omega _s^r\tau +\phi _e\big)}+ {\mathcal P}_a^a\cos {\big(\omega _s^a\tau +\phi _a\big)}, \\[-12pt]\nonumber \end{align}

(3.6e) \begin{align} \overline {\varpi }_{2s}(\tau )&=\dfrac {\varpi _{2s}}{\varepsilon _r k_y c_2}= {\mathcal O}_r^r\cos {\big(\omega _s^r\tau +\phi _r\big)+ \mathcal O}_r^a\cos {\big(\omega _s^a\tau +\phi _a\big)}. \end{align}

The postshock mode amplitudes ${\mathcal U}_{r,a}^{r,a}$ , ${\mathcal V}_{r,a}^{r,a}$ , ${\mathcal D}_{e,a}^{r,a}$ , ${\mathcal P}_{a}^{r,a}$ and ${\mathcal O}_{r}^{r,a}$ are unknowns functions of $\theta ^r$ and $\theta ^a$ to be determined. Note that in (3.5) and (3.6), the transverse coordinate dependence is omitted, $\sim\! \cos (k_y y)$ for the longitudinal and $\sim\! \sin (k_y y)$ for the transverse velocity field, as it is not affected by the shock passage due to the periodic symmetry in the transverse direction. To distinguish between the different contributions, the subscripts $r,e,a$ refer to the nature of the downstream perturbation, rotational, entropic or acoustic, respectively, while the superscripts $r,a$ refer to the nature of the upstream mode generating the corresponding perturbation, rotational entropic and acoustic, respectively. As such, the rotational streamwise amplitude ${\mathcal U}_{r}(\theta ^r,\theta ^a)$ is conveniently decomposed into components associated with the upstream rotational-entropic wave ${\mathcal U}_{r}^r={\mathcal U}_{r}^r(\theta ^r)$ and the acoustic mode ${\mathcal U}_{r}^a= {\mathcal U}_{r}^a(\theta ^a)$ , with the superposition of modes rendering ${\mathcal U}_{r}={\mathcal U}_{r}^r+{\mathcal U}_{r}^a$ . Likewise, ${\mathcal U}_{a}(\theta ^r,\theta ^a)={\mathcal U}_{a}^r(\theta ^r)+{\mathcal U}_{a}^a(\theta ^a)$ , ${\mathcal D}_{e}(\theta ^r,\theta ^a)={\mathcal D}_{e}^r(\theta ^r)+{\mathcal D}_{e}^a(\theta ^a)$ , and so on.

In the compressed gas reference frame, $x_c=x-(u_1-u_2)t$ , the rotational-entropic disturbances remain frozen (steady functions) to the fluid particles, while acoustic perturbations propagate governed by the corresponding sound wave equation

(3.7) \begin{equation} \dfrac {\partial ^2\overline {p}_2}{\partial \tau ^2}=\dfrac {\partial ^2\overline {p}_2}{\partial \overline {x}_c^2}-\overline {p}, \end{equation}

where $\overline {x}_c=k_y x_c$ , thereby behaving similar to upstream disturbances described in (3.2d ). To compute the long-time amplitudes of the different modes, the perturbed RH equations must be employed, along with the linearised Euler equations. An additional condition associated with the shock supporting mechanism is needed; in this case, as usual, the isolated shock boundary condition is employed far downstream of the shock. Therefore, all disturbances are generated at the shock front and conveyed (rotational entropic) or radiated (acoustic) downstream.

The boundary condition at the shock front is derived from the linearised RH jump conditions assuming that (a) the thickness of the vibrational non-equilibrium region $\ell _T$ is much smaller than the inverse of the transverse wavenumber $k_y^{-1}$ and (b) the displacement of the shock $\xi _s=\xi _s(y,t)$ from its mean, flat shape is much smaller than $k_y^{-1}$ . In these limits, at any transverse coordinate $\overline {y}=k_y y$ , the RH jump conditions can be applied at the mean shock front location $\overline {x}_c= {\mathcal M}_2\tau$ (measured from the compressed gas reference frame), and can be linearised about the mean thermochemical equilibrium postshock gas state $\mathcal{R} = \rho _2/\rho _1$ , $\mathcal{T} = T_2/T_1$ and $\mathcal{M}_2$ , along with the normalised RH-slope parameters $\varGamma$ , $\varGamma _\rho$ and $\varGamma _{\!p}$ , calculated in § 2, thereby yielding

(3.8a) \begin{align} \frac {\textrm {d} \overline {\xi }_s}{\textrm {d} \tau } &= \frac { {\mathcal R} \left (1-\varGamma \right )}{2 {\mathcal M}_2 \left ({\mathcal R}-1\right )}\overline {p}_{2s} - \frac {\mathcal{M}_2\mathcal{R}^2\left (1- 2\mathcal{R}^{-1}-\varGamma _\rho \right )}{2\left (\mathcal{R} - 1\right )} \overline {\rho }_{1s} + \dfrac {1}{\beta }\overline {u}_{1s} - \dfrac {1 - \varGamma _{\!p}}{2 \beta ^2 \mathcal{M}_2 (\mathcal{R} - 1)} \overline {p}_{1s}, \\[-12pt]\nonumber \end{align}

(3.8b) \begin{align} \overline {u}_{2s} &=\frac {1+\varGamma }{2 {\mathcal M}_2}\overline {p}_{2s} - \frac {\mathcal{M}_2\mathcal{R}\left (1 +\varGamma _\rho \right )}{2} \overline {\rho }_{1s} + \dfrac {1}{\beta }\overline {u}_{1s} - \dfrac {1 + \varGamma _{\!p}}{\beta ^2\mathcal{M}_2\mathcal{R}} \overline {p}_{1s}, \\[-12pt]\nonumber \end{align}

(3.8c) \begin{align} \overline {v}_{2s} &= \dfrac {1}{\beta } \overline {v}_{1s}-{\mathcal M}_2 \left ({\mathcal R}-1\right )\frac {\partial \overline {\xi }_s}{\partial \overline {y}}, \\[-12pt]\nonumber \end{align}

(3.8d) \begin{align} \overline {\rho }_{2s} &= \frac {\varGamma }{{\mathcal M}_2^2}\overline {p}_{2s} - \mathcal{R}\varGamma _\rho \overline {\rho }_{1s} - \dfrac {\varGamma _{\!p}}{\beta ^2 \mathcal{M}_2^2 \mathcal{R}} \overline {p}_{1s}. \end{align}

In the above equations, $\overline {\xi }_s = k_y\xi _s/\varepsilon _r$ is the dimensionless shock displacement, and $\beta = c_2/c_1$ is the ratio of the sound speeds across the shock. Additionally, to ensure that thermochemical equilibrium is reached, the thickness of the vibrational non-equilibrium region $\ell _T$ must be much smaller than the inverse of the transverse wavenumber $k_y^{-1}$ , i.e. $k_y \ell _T \ll 1$ . If upstream density deviations are positively correlated with rotational disturbances, i.e. $\varepsilon _e/\varepsilon _r\gt 0$ , they contribute to reduce the shock corrugation amplitude and postshock streamwise velocity perturbation. The opposite occurs for negatively correlated disturbances. Note also that, for strong shocks ${\mathcal M}_1\gg 1$ , density disturbances (either entropic or acoustic) become more relevant because the influence of velocity and pressure perturbations diminishes as they are proportional to $\beta ^{-1}\ll 1$ and $\beta ^{-2}\ll 1$ , respectively. These trends are captured with opposite signs compared with studies such as those of Mahesh et al. (Reference Mahesh, Lele and Moin1997), Sinha (Reference Sinha2012) and Quadros et al. (Reference Quadros, Sinha and Larsson2016a , Reference Quadros, Sinha and Larssonb ), owing to the orientation of our coordinate system: in the present work, the $x$ -axis is directed from downstream to upstream, opposite to the flow direction (see figure 1). As a result, a positive correlation $\varepsilon _e/\varepsilon _r \gt 0$ in their formulation corresponds to $\varepsilon _e/\varepsilon _r \lt 0$ in ours.

In the long-time regime, the factor $\partial \bar {\xi }_s/\partial \tau$ drops two frequency factors $\omega _s^r$ and $\omega _s^a$ corresponding to the different sources of the forcing terms: rotational entropic and acoustic. To obtain the value of the perturbation amplitudes downstream, linear superposition can be applied to independently solve (3.7)–(3.8) for the two upstream sources of disturbances.

3.2. Amplification of turbulent kinetic energy, turbulent intensity and turbulent Reynolds number across the shock

The weak isotropic turbulence in the preshock gas can be modelled as a linear superposition of incident rotational-entropic and acoustic waves, with the corresponding amplitudes $\varepsilon _r$ , $\varepsilon _e$ and $\varepsilon _a$ varying according to an isotropic energy spectrum. In three-dimensional turbulence, the spectrum is related to the modal amplitude by $E(k) = 4\pi k^2 \varepsilon ^2(k)$ , whereas in two-dimensional turbulence it is given by $E(k) = 2\pi k \varepsilon ^2(k)$ . By assuming isotropy – which implies that the probability of an incident wave having a specific orientation is proportional to the solid angle in three-dimensional space – the root mean square of the velocity and vorticity fluctuations in the preshock gas can be determined. However, to compute the statistical averages downstream, the correlation between the different modes must be specified.

In what concerns the relation between entropic and rotational disturbances, Mahesh et al. (Reference Mahesh, Lele and Moin1997) has addressed the two cases $\varepsilon _e\sim \varepsilon _r /{\mathcal M}_1$ and $\varepsilon _e \sim -\varepsilon _r /{\mathcal M}_1$ with the same spectral density distribution for each type of fluctuations. If $\varepsilon _e$ had been chosen to be proportional to $\varepsilon _r {\mathcal M}_1$ instead of $\varepsilon _r /{\mathcal M}_1$ , as it would be for the disturbances suggested by Morkovin (Reference Morkovin1962), then the entropic contribution would be much more important at high Mach numbers. In turbulent boundary layers, where a preferential direction is imposed by the free-stream flow, it is natural to relate the amplitude of entropic disturbances to that of the streamwise velocity perturbations, as previously noted. In contrast, for cases like those considered here – where the turbulence is homogeneous and no preferred direction exists – entropic disturbances should be compared with the amplitude of the full velocity perturbation vector, rather than only its longitudinal component. A related parameter characterises the correlation between entropic and solenoidal disturbances, which can be generalised in terms of the relative phase in a normal mode analysis, i.e. $\varepsilon _e \sim \varepsilon _r \exp (i\,\phi _{er})$ . The choice of this phase $\phi _{er}$ influences the final statistical averages, as it determines the degree of constructive or destructive interference between solenoidal turbulence amplification and baroclinic turbulence generation across the shock.

Similarly to the works of Mahesh et al. (Reference Mahesh, Lele and Moin1997), Veera & Sinha (Reference Veera and Sinha2009) and Quadros et al. (Reference Quadros, Sinha and Larsson2016a , Reference Quadros, Sinha and Larssonb ), we define the ratio between entropic and vortical fluctuations as $\chi \exp (i\ \phi _{re}) = \varepsilon _e/\varepsilon _r$ , whose averaged value

(3.9) \begin{equation} \langle \chi \rangle = \dfrac {\big\langle \delta \rho _1^e \delta u_1^r\big\rangle }{\big\langle \delta u_1^r \delta u_1^r \big\rangle } \dfrac {\langle c_1 \rangle }{\langle \rho _1 \rangle } \end{equation}

can be extracted from the numerical characterisation of the upstream turbulence flow.

As for the acoustic influence, the correlation with vorticity disturbances differs due to their distinct nature as travelling waves, resulting in no net value from the product of these two contributions. Therefore, the impact of the acoustic waves in the TKE can be evaluated separately, with the total contribution being simply computed as the weighted sum of the respective parts. For this, we define the dimensionless variable $\eta =\varepsilon _a^2/\varepsilon _r^2$ , which, in terms of turbulent variables, is computed as

(3.10) \begin{equation} \eta = \dfrac {\text{TKE}_{1a}}{\text{TKE}_{1r}}=\dfrac { \big\langle \delta u_{1a}^2 \big\rangle + \big\langle \delta v_{1a}^2 \big\rangle + \big\langle \delta w_{1a}^2 \big\rangle }{\big\langle \delta u_{1r}^2\big\rangle +\big\langle \delta v_{1r}^2 \big\rangle +\big\langle \delta w_{1r}^2 \big\rangle }, \end{equation}

where $\langle {\delta u}_{1a}^2 \rangle , \langle {\delta v}_{1a}^2 \rangle$ and $\langle {\delta w}_{1a}^2 \rangle$ denote the averaged dimensional value of the three orthogonal contributions of the acoustic velocity field ahead of the shock. Correspondingly, $\langle {\delta u}_{1r}^2 \rangle , \langle {\delta v}_{1r}^2 \rangle$ and $\langle {\delta w}_{1r}^2 \rangle$ refer to upstream rotational perturbations.

The analysis follows a similar approach to that outlined in Wouchuk et al. (Reference Wouchuk, Huete Ruiz de Lira and Velikovich2009), Huete et al. (Reference Huete, Velikovich and Wouchuk2011, Reference Huete, Wouchuk and Velikovich2012b , Reference Huete, Cuadra, Vera and Urzay2021), Cuadra et al. (Reference Cuadra, Vera, Di Renzo and Huete2023) and Cuadra (Reference Cuadra2023). For the sake of conciseness, only the main steps will be presented here, without repeating the detailed procedures from those references. The isotropy condition implies that the probability of the incident wave having a specific orientation angle is proportional to the solid angle, given by $\sin \theta\, \textrm {d}\theta\, \textrm {d}\varphi /(4\pi )$ . Here, $\varphi$ is defined as part of the extension from a two-dimensional velocity field to a three-dimensional configuration, as explained in Huete et al. (Reference Huete, Cuadra, Vera and Urzay2021).

The total TKE amplification factor, $K$ , is defined as follows:

(3.11) \begin{equation} K=\dfrac {\text{TKE}_{2}}{\text{TKE}_{1}}=\dfrac { \big\langle \delta u_{2}^2 \big\rangle + \big\langle \delta v_{2}^2 \big\rangle + \big\langle \delta w_{2}^2 \big\rangle }{\big\langle \delta u_{1}^2 \big\rangle +\big\langle \delta v_{1}^2\big\rangle + \big\langle \delta w_{1}^2 \big\rangle }=\beta ^2\dfrac { \big\langle \overline {u}_{2}^2 \big\rangle + \big\langle \overline {v}_{2}^2 \big\rangle + \big\langle \overline {w}_{2}^2 \rangle }{\big\langle \overline {u}_{1}^2 \big\rangle + \big\langle \overline {v}_{1}^2 \big\rangle + \big\langle \overline {w}_{1}^2 \big\rangle }, \end{equation}

where the averaged dimensional values of the perturbation kinetic energy associated with the velocity components behind and ahead of the shock are denoted with subscripts $2$ and $1$ , respectively. The factor $\beta ^2$ appears because dimensionless velocity perturbations are scaled with the local speed of sound, which change across the shock.

Let us start with the computation of the TKE generated by rotational-entropic disturbances $K^r$ . This factor can be conveniently split into its longitudinal ( $K_L^r$ ) and transverse ( $K_T^r$ ) contributions through

(3.12) \begin{equation} K^r=\dfrac {1}{3}\left (K_L^r+2 K_T^r\right )=\dfrac {\beta ^2}{3}\left (L^r+2 T^r\right ), \end{equation}

in which either the longitudinal $L^r=K_L^r/\beta ^2$ and the transverse $T^r=K_T^r/\beta ^2$ contributions to the TKE have been normalised with the ratio of speed of sound squared $\beta ^2$ . In addition, they can be separated into their acoustic and rotational contributions downstream, namely $L^r=L_r^r+L_a^r$ and $T^r=T_r^r+T_a^r$ . With use made of the isotropy condition, and following the steps detailed in Huete et al. (Reference Huete, Cuadra, Vera and Urzay2021), each of these contributions is written as an explicit analytical integral equation. For example, the functions

(3.13a) \begin{align} L_r^r &= \dfrac {3}{2}\int _0^{\pi /2}\big [ \mathcal{U}_r^{r}(\theta ^{r})\big ]^2\sin ^{3}\theta ^{r} \textrm {d}\theta ^r, \\[-12pt]\nonumber \end{align}

(3.13b) \begin{align} L_a^r &= \dfrac {3}{2}\int _0^{\theta _{r}^*}\big [ \mathcal{U}_a^{r}(\theta ^{r})\big ]^2\sin ^{3}\theta ^{r} \textrm {d}\theta ^r , \end{align}

are normalised so that $L_r^r=L_a^r=1$ for ${\mathcal U}_r^r={\mathcal U}_a^r=1$ . In writing (3.13) the symmetry condition $\int _{0}^{\pi /2} \text{fun}(\theta ^r)\, \textrm {d}\theta ^r= 2\int _0^{\pi /2} \text{fun}(\theta ^r)\, \textrm {d}\theta ^r$ for the rotational contribution and $\int _{-\theta _{r}^*}^{\theta _{r}^*} \text{fun}(\theta ^r)\, \textrm {d}\theta ^r= 2\int _0^{\theta _{r}^*} \text{fun}(\theta ^r)\, \textrm {d}\theta ^r$ for the acoustic part have been used to reduce the integration domain. The angle given by

(3.14) \begin{equation} \tan \theta ^{r*}=\dfrac {{\mathcal M}_2{\mathcal R}}{\sqrt {1-{\mathcal M}_2^2}} \end{equation}

distinguishes the short wavelength (high-frequency) and long wavelength (low-frequency) regimes, corresponding to $\theta ^r\lt \theta ^{r*}$ and $\theta ^r\gt \theta ^{r*}$ , respectively. For the latter, the acoustic disturbances generated behind the shock are evanescent waves, with the amplitude decaying exponentially from the distance to the shock. For this reason, the acoustic contribution to the longitudinal kinetic energy in the asymptotic far field only includes the contribution of the short-wavelength regime.

The dependence of the asymptotic amplitudes on the upstream rotational wave-vector angle, ${\mathcal U}_r^r(\theta ^r)$ and ${\mathcal U}_a^r(\theta ^r)$ , is obtained through a standard normal mode analysis. This procedure – namely, the decomposition of small perturbations into harmonic, small-amplitude waves – follows the methodology presented in Huete et al. (Reference Huete, Cuadra, Vera and Urzay2021), which itself builds upon the seminal work of Ribner (Reference Ribner1954a ). In the present formulation, the terms associated with the perturbation of the RH curve are extended to account for the variations arising under hypersonic conditions. The resulting amplitudes, which also depend on the coupling between rotational and entropic disturbances through the factor $\chi$ , in a form analogous to that obtained in Huete, Sánchez & Williams (Reference Huete, Sánchez and Williams2014) for exothermic reacting shocks, are expressed as a piecewise function separated by the critical angle $\theta ^{r*}$ . To avoid repeating this standard analysis, the corresponding expressions are provided in the Combustion Toolbox (Cuadra et al. Reference Cuadra, Huete and Vera2025b ).

The integral expressions for the transverse contributions generated by the upstream rotational-entropic disturbances are given by

(3.15a) \begin{align} T_r^r&= \dfrac {3}{4}\int _0^{\pi /2}\left [{\mathcal V}_r^r(\theta ^r)\right ]^2\sin ^3\theta ^r \textrm {d}\theta ^r + \dfrac {3}{4}\dfrac {1}{\beta ^2}, \\[-12pt]\nonumber \end{align}

(3.15b) \begin{align} T_a^r&= \dfrac {3}{4}\int _0^{\theta _{r}^*}\left [{\mathcal V}_a^r(\theta ^r)\right ]^2\sin ^3\theta ^r \textrm {d}\theta ^r, \end{align}

where the asymptotic amplitudes, ${\mathcal V}_r^r(\theta ^r)$ and ${\mathcal V}_a^r(\theta ^r)$ are obtained using the same normal mode analysis employed for the corresponding longitudinal components.

The total TKE defined in (3.11) includes the part generated by upstream rotational-entropic disturbances ( $K^r$ ), presented in (3.12), and the TKE generated by the acoustic fluctuations ahead of the shock ( $K^a$ ). The latter, following the analysis detailed in Huete et al. (Reference Huete, Wouchuk and Velikovich2012b ), may be defined in similar terms

(3.16) \begin{equation} K^a=\dfrac {1}{3}\left (K_L^a+2 K_T^a\right )=\dfrac {\beta ^2}{3}\left (L^a+2 T^a\right ), \end{equation}

that is also split into its longitudinal ( $K_L^a$ and $L^a$ ) and transverse ( $K_T^a$ and $T^a$ ) contributions. These functions require special consideration because the upstream acoustic field is not symmetric about the $\theta =0$ axis. This asymmetry arises from the nature of acoustic waves as travelling waves, distinguishing those moving in the direction of shock propagation from those directed toward the shock surface. In what concerns the longitudinal kinetic energy generated by the upstream acoustic field, a similar decomposition is employed $L^a=L_r^a+L_a^a$ , with each part being explicitly written as

(3.17a) \begin{align} L_r^a&= \dfrac {3}{4}\int _{0}^{\theta ^{a*}}\left [{\mathcal U}_r^a(\theta ^a)\right ]^2\sin ^3\theta ^a \textrm {d}\theta ^a + \dfrac {3}{4}\int _{\theta ^{a\dagger }}^{\pi }\left [{\mathcal U}_r^a(\theta ^a)\right ]^2\sin ^3\theta ^a \textrm {d}\theta ^a,\\[-12pt]\nonumber \end{align}

(3.17b) \begin{align} L_a^a&= \dfrac {3}{4}\int _{\theta ^{a*}}^{\theta ^{a\dagger }}\left [{\mathcal U}_a^a(\theta ^a)\right ]^2\sin ^3\theta ^a \textrm {d}\theta ^a. \\[9pt] \nonumber \end{align}

In this case, the critical acoustic wavenumber angles is readily obtained by the implicit relationship

(3.18) \begin{equation} \tan \theta ^a=\left (1-\dfrac {1}{{\mathcal M}_1\cos \theta ^a}\right )\dfrac {{\mathcal M}_2{\mathcal R}}{\sqrt {1-{\mathcal M}_2^2}}, \end{equation}

with $\theta ^{a*}$ and $\theta ^{a\dagger }$ being used for the positive and negative propagation directions, respectively, defined in the domain $0\leqslant \theta ^a \leqslant \pi$ ; different values of the forward moving or backwards moving acoustic wave produce similar shock oscillation frequencies depending on the shock propagation speed ${\mathcal M}_1$ . In this case, the range $\theta ^{a*}\leqslant \theta ^a\leqslant \theta ^{a\dagger }$ corresponds to the short-wavelength regime, while the domain bounded by $0\leqslant \theta ^a\leqslant \theta ^{a*}$ and $\theta ^{a\dagger }\leqslant \theta ^a\leqslant \pi$ refer to the long-wavelength regime, referring to evanescent or non-decaying acoustic waves behind the shock, respectively. In the limit ${\mathcal M}_1\gg 1$ , this effect may become negligible as the acoustic wave speed is much smaller than the shock speed.

The transverse component of the TKE generated by turbulent acoustic perturbations, $T^a=T_r^a+T_a^a$ , is expressed in integral expressions for each term

(3.19a) \begin{align} T_r^a&= \dfrac {3}{8}\int _{0}^{\theta ^{a*}}\left [{\mathcal V}_r^a(\theta ^a)\right ]^2\sin ^3\theta ^a \textrm {d}\theta ^a + \dfrac {3}{8}\int _{\theta ^{a\dagger }}^{\pi }\left [{\mathcal V}_r^a(\theta ^a)\right ]^2\sin ^3\theta ^a \textrm {d}\theta ^a, \\[-12pt]\nonumber \end{align}

(3.19b) \begin{align} T_a^a&= \dfrac {3}{8}\int _{\theta ^{a*}}^{\theta ^{a\dagger }}\left [{\mathcal V}_a^a(\theta ^a)\right ]^2\sin ^3\theta ^a \textrm {d}\theta ^a, \end{align}

for the rotational and acoustic contributions downstream, respectively. Here, the dependence of the amplitudes on the upstream acoustic wave-vector angle, ${\mathcal U}_r^a(\theta ^a)$ and ${\mathcal U}_a^a(\theta ^a)$ for the longitudinal contribution, and ${\mathcal V}_r^a(\theta ^a)$ and ${\mathcal V}_a^a(\theta ^a)$ for the transverse components, can be obtained using the same normal mode analysis described in Huete et al. (Reference Huete, Wouchuk and Velikovich2012b ), and is conveniently implemented in the Combustion Toolbox (Cuadra et al. Reference Cuadra, Huete and Vera2025b ) for reference.

With the partial contributions defined, the total kinetic energy amplification factor (3.11), and its partial contributions, can be written as a linear superposition of the upstream rotational and acoustic contributions, namely

(3.20) \begin{equation} K=\dfrac {K^r+ \eta K^a}{1+\eta },\quad K_L=\beta ^2\dfrac {L^r+ \eta L^a}{1+\eta },\quad K_T=\beta ^2\dfrac {T^r+ \eta T^a}{1+\eta }. \end{equation}

Therefore, the relation between the three modes characterising the upstream turbulent perturbations, defined as $\varepsilon _e/\varepsilon _r$ and $\varepsilon _a/\varepsilon _r$ in the mono-frequency configuration, are included in the total TKE amplification factor $K$ through the coefficients $\chi$ , which is internal to the computation of $K^r$ , and $\eta$ , which is factored out in the weighted sum of contribution for $K^a$ , $L^a$ and $T^a$ .

Figure 3. Far-field TKE amplification ratio $K$ (a) and its acoustic contribution $K_a$ (b) predicted by LIA as a function of the preshock Mach number $\mathcal{M}_1$ for different intensities of dilatational-to-solenoidal TKE $\eta$ . The calculations are provided for the thermally perfect (solid lines) and calorically perfect (dashed lines) models.

Figure 3(a) presents the TKE amplification ratio, $K$ , as a function of the preshock Mach number, $\mathcal{M}_1$ , for various intensities of dilatational-to-solenoidal TKE, $\eta$ . The results are shown for both the thermally perfect (solid lines) and calorically perfect (dashed lines) models. It is observed that dilatational disturbances make a significant contribution to the turbulence amplification factor for sufficiently strong shock intensities, $\mathcal{M}_1\geqslant 5$ , owing to the dominant mechanisms of vorticity tilting and compressive amplification. To understand this behaviour, it is necessary to focus on the amplification of kinetic energy due solely to upstream acoustic energy, denoted by $K^a$ . This is illustrated in figure 6(b) of Huete et al. (Reference Huete, Wouchuk and Velikovich2012b ), where the solid line crosses unity at two distinct Mach numbers, ${\mathcal M}_1=1.11$ and ${\mathcal M}_1=2.35$ . Under these conditions, there is no net change in TKE due to upstream acoustic disturbances, irrespective of the amplitude $\eta$ , as indicated by the intersection of the different curves in figure 3(a) at these Mach numbers. In the hypersonic regime, the inclusion of compressibility effects in the upstream turbulence results in an increase in TKE, with an asymptotic contribution that grows with ${\mathcal M}_1^2$ . Figure 3(b) presents the downstream TKE associated with acoustic disturbances $K_a$ as a function of $\mathcal{M}_1$ . It is evident that the downstream acoustic energy increases when dilatational disturbances are present in the upstream flow, although their relative contribution remains relatively modest, with $K_a/K\lesssim 5\,\%$ .

Figure 4. Far-field streamwise $K_L$ (a) and transverse $K_T$ (b) components of the TKE amplification ratio predicted by LIA as a function of the preshock Mach number $\mathcal{M}_1$ for different intensities of dilatational-to-solenoidal TKE $\eta$ . The calculations are provided for the thermally perfect (solid lines) and calorically perfect (dashed lines) models.

When vibrational effects are taken into account, a similar trend is observed, albeit with slightly different Mach numbers at the intersection points. For the case ${\mathcal M}_1 = 5$ , the effect of considering a thermally versus calorically perfect gas is not qualitatively significant. However, this does not hold when separating the streamwise and transverse contributions, denoted as $K_L = \langle \delta u_2^2 \rangle / \langle \delta u_1^2 \rangle$ and $K_T = \langle \delta v_2^2 + \delta w_2^2 \rangle / \langle \delta v_1^2 + \delta w_1^2 \rangle$ , respectively. These are shown in figures 4(a) and 4(b) as functions of the Mach number $\mathcal{M}_1$ . For the streamwise component of the kinetic energy amplification factor, the consideration of a temperature-dependent $\gamma (T)$ leads to a reduction in $K_L$ , whereas the opposite effect is observed for $K_T$ . These opposing trends largely offset each other when evaluating the total TKE computed in figure 3. The reason for the increase of the transverse contribution is analysed in detail in Huete et al. (Reference Huete, Cuadra, Vera and Urzay2021); it is associated with the amplification of the density compression ratio that ultimately leads to an increase in the local flow deflection of the flow.

3.3. One-dimensional spectral formulation

This subsection presents the theoretical framework used to compute the one-dimensional spectra of various turbulent quantities using linear theory. The methodology provides a mode-by-mode quantification of shock-induced spectral amplification given an upstream turbulence spectrum, similar to the approach outlined by Ribner (Reference Ribner1986, Reference Ribner1987) for shock waves and Jackson, Hussaini & Ribner (Reference Jackson, Hussaini and Ribner1993) and Cuadra, Huete & Vera (Reference Cuadra, Huete and Vera2020) for detonation waves. The present formulation builds upon this foundation by expressing the downstream spectrum as a weighted composition of solenoidal and dilatational contributions, each modulated by their respective LIA amplification factors, consistent with the decomposition introduced in the previous subsections. The resulting spectra are later evaluated and directly compared with DNS data in § 5.4.

To gain further insight into the structure of the turbulent field that develops far downstream of the shock, the one-dimensional power spectra of the relevant fluctuating quantities can be evaluated through the integral functions

(3.21) \begin{equation} \varPhi _{j}^r(k^r) = N^r\int _0^{\pi /2} {\mathcal G}_{\!j}^r(\theta ^r)\varepsilon _r^2(k^r) [k^r]^2 \sin ^3\theta ^r\ \mathrm{d}\theta ^r, \end{equation}

where ${\mathcal G}_{\!j}^r(\theta ^r)$ is an auxiliary function that represents the corresponding single-frequency perturbation amplitudes behind the shock produced by the rotational disturbances upstream, namely $ [\mathcal{P}_a^r(\theta ^r) ]^2$ , $ [\mathcal{U}_r^r(\theta ^r) ]^2 + [\mathcal{U}_a^r(\theta ^r) ]^2$ , $ [\mathcal{V}_r^r(\theta ^r) ]^2 + [\mathcal{V}_a^r(\theta ^r) ]^2$ and $ [\mathcal{D}_e^r(\theta ^r) ]^2 + [\mathcal{D}_a^r(\theta ^r) ]^2$ for $j=p,u,v$ and $\rho$ , respectively. The subscript $j=o$ will refer to preshock conditions, and therefore ${\mathcal G}_o^r=1$ for that case. Those associated with the acoustic perturbation upstream can be equally defined

(3.22) \begin{equation} \varPhi _{\!j}^a(k^a) = N^a \int _{0}^{\pi } {\mathcal G}_{\!j}^a(\theta ^a) \varepsilon _a^2(k^a) [k^a]^2\sin ^3\theta ^a\ \mathrm{d}\theta ^a, \end{equation}

where the auxiliary function ${\mathcal G}_{\!j}^a(\theta ^a)$ now reads as $ [\mathcal{P}_a^a(\theta ^a) ]^2$ , $ [\mathcal{U}_r^a(\theta ^a) ]^2+ [\mathcal{U}_a^a(\theta ^a) ]^2$ , $ [\mathcal{V}_r^a(\theta ^a) ]^2 + [\mathcal{V}_a^a(\theta ^a) ]^2$ and $ [\mathcal{D}_e^a(\theta ^a) ]^2+ [\mathcal{D}_a^a(\theta ^a) ]^2$ for $j=p,u,v$ and $\rho$ , respectively. Likewise, $j=o$ refers to the upstream spectrum, ${\mathcal G}_o^a=1$ .

The total weighted contribution reads as

(3.23) \begin{equation} \varPhi _{\!j}(k^r, k^a) = \dfrac {\varPhi _{\!j}^r(k^r) + \eta \, \varPhi _{\!j}^a(k^a)}{1 + \eta },\quad j = p, u, v, \rho , o, \end{equation}

where the normalisation parameters $N^r$ and $N^a$ are defined to ensure unit energy content

(3.24a) \begin{align} [N^r]^{-1} &= \int _0^{\infty }\int _0^{\pi /2}\varepsilon _r^2(k^r) [k^r]^2\sin ^3\theta ^r\ \mathrm{d} \theta ^r \mathrm{d} k^r, \\[-12pt]\nonumber \end{align}

(3.24b) \begin{align} [N^a]^{-1} &= \int _0^{\infty }\int _0^{\pi }\varepsilon _a^2(k^a) [k^a]^2\cos ^2\theta ^a\sin \theta ^a\ \mathrm{d} \theta ^a \mathrm{d} k^a. \\[9pt] \nonumber \end{align}

Note that in (3.21)–(3.22), the integral domain in the polar angle coordinates $\theta ^r$ and $\theta ^a$ , may not necessarily cover $[0,\pi /2]$ and $[0,\pi ]$ , as acoustic disturbances only contribute in the far field when $\theta ^r\leqslant \theta ^{r*}$ and/or $\theta ^{a*}\leqslant \theta ^a\leqslant \theta ^{a\dagger }$ , respectively. This has been omitted for the sake of simplicity in presenting the integral spectrum equations.

The spectral amplitudes of vortical-entropic and acoustic fluctuations are related to their respective three-dimensional isotropic spectral densities, namely

(3.25a) \begin{align} \varPhi _{11}^r(\boldsymbol{k^r}) &= \dfrac {E^r(k^r) \sin ^2\theta ^r}{4\pi [k^r]^2} = \varepsilon _r^2(k^r) \sin ^2\theta ^r, \\[-12pt]\nonumber \end{align}

(3.25b) \begin{align} \varPhi _{11}^a(\boldsymbol{k^a}) &= \dfrac {E^a(k^a)\cos ^2\theta ^a}{8\pi [k^a]^2} = \dfrac {1}{2}\varepsilon _a^2(k^a) [k^a]^2 \cos ^2\theta ^a, \end{align}

where $E^r(k^r)$ and $E^a(k^a)$ denote the isotropic energy spectrum of the respective component. The former can be prescribed analytically, for example, using the von Kármán (VK) model, which captures the Kolmogorov scaling in the inertial subrange ( $\sim\! k^{-5/3}$ ) as well as the low-wavenumber behaviour ( $\sim\! k^4$ )

(3.26) \begin{equation} E^r(k^r) \sim \dfrac {\left (k^r/k_0\right )^4}{\left [1 + \left (k^r/k_0\right )^2\right ]^{17/6}}, \end{equation}

where $k_0$ represents the scale frequency parameter defining the energy-containing subrange. In contrast, there is no equivalent closed-form expression for the acoustic contribution, $E^a(k^a)$ , as it typically arises as a correction to the rotational spectrum and cannot, in general, be represented independently. However, the VK expression will be used only as a reference, since it does not accurately reproduce the energy distribution of small-scale turbulent fluctuations, which follow a steeper decay power law. Alternatively, both spectra can be extracted from DNSs by applying a Helmholtz decomposition to the velocity field, as described in Petersen & Livescu (Reference Petersen and Livescu2010). Note that in this case, the solenoidal and dilatational content are already prescribed by the turbulence generator. Therefore, the kinetic energy spectrum is simply given by the sum of the solenoidal and dilatational contributions.

4.1. Conservation equations and numerical method

The DNSs have been carried out with the Hypersonics Task-based Research (HTR) solver (Di Renzo, Fu & Urzay Reference Di Renzo, Fu and Urzay2020; Di Renzo & Pirozzoli Reference Di Renzo and Pirozzoli2021; Di Renzo Reference Di Renzo2022), numerically integrating the single-component compressible Navier–Stokes equations. In this formulation, the temperature-dependent viscosity is evaluated with Sutherland’s law specialised for air, whereas the thermal conductivity of the gas is determined using a constant Prandtl number of 0.71. The frozen chemical composition of the gas comprises 79 % molecular nitrogen and 21 % molecular oxygen on a molar basis. Two thermodynamic models are considered in these calculations, pertaining to the description of the internal-energy modes of air. The first model treats air as a calorically perfect diatomic gas whose constant-pressure heat capacity is constant and equal to $c_{\!p} = 7/2 \mathcal{R}^0$ , whereas the second represents air as a thermally perfect gas with equilibrium vibrational excitation. For the latter model, the constant-pressure heat capacity is evaluated using a nine-coefficient polynomial particularised for dry air absent argon (McBride Reference McBride2002). These two models represent the two ends of a broad range of regimes for the interaction of thermal relaxation phenomena with STI. The first (second) model is applicable when the thermal relaxation is infinitely slower (faster) than turbulence and convective time scales. Intermediate scenarios, where the thermal relaxation of the gas molecules is a finite-rate process compared with hydrodynamics, will likely take place in realistic conditions. However, such a mixed scenario cannot be directly predicted by LIA, and therefore, for the sake of having the most fare comparison across linear and nonlinear descriptions of the hypersonic STI, we have decided to execute the DNS only in these two limits and to defer a quantification of the impact of thermodynamic non-equilibrium to future works.

The conservation equations corresponding to mass, momentum and total energy are discretely solved using conservative finite differences on Cartesian grids. The inviscid fluxes are discretised using a hybrid reconstruction scheme that utilises a skew-symmetric formulation (Pirozzoli Reference Pirozzoli2010) in smooth regions of the flow and the (Jiang & Shu Reference Jiang and Shu1996) fifth-order weighted essentially non-oscillatory (WENO5-JS) reconstruction scheme in the vicinity of discontinuities. The switch between the two schemes is governed by a directional shock sensor based on the targeted essentially non-oscillatory (TENO) (Fu, Hu & Adams Reference Fu, Hu and Adams2016) stencil selection process (Williams, Di Renzo & Moin Reference Williams, Di Renzo and Moin2022). The resulting scheme is well suited for numerical solution of the STI problem (Johnsen et al. Reference Johnsen2010), as the hybridised method provides very low numerical dissipation away from discontinuities due to the conservation properties of the selected skew-symmetric numerical scheme, combined with a sharp and stable representation of shock waves with the high-resolution shock-capturing method. The diffusion fluxes are discretised using a second-order central formulation written in conservative form. Based on the semi-discretisation of the conservation equations with the aforementioned spatial fluxes, the resulting set of ordinary differential equations is advanced in time using a third-order strong stability-preserving Runge–Kutta scheme (Gottlieb, Shu & Tadmor Reference Gottlieb, Shu and Tadmor2001).

Figure 5. Schematic of the computational set-up showcasing three-dimensional isotropic turbulence generated by a triply periodic box and advecting at a mean Mach number $\mathcal{M}_1$ toward the STI domain.

The computational framework for the DNS consists of two distinct calculations executed concurrently as depicted in figure 5. The first simulation, functioning as a turbulence generator, comprises a triply periodic cubic computational domain with dimensions $2 \pi \times 2\pi \times 2\pi$ wherein HIT is forced to sustain TKE and has the mean convection velocity to be imposed in the preshock conditions. The HIT domain is discretised with $800^3$ uniformly distributed points. The latter calculation corresponds to the STI itself, consisting of a computational domain with dimensions $4 \pi \times 2\pi \times 2\pi$ . The computational domain is periodic in the transverse directions, namely $y$ and $z$ , and it is discretised with $1600\,\times \,800^2$ points in the streamwise and transverse directions, respectively. The grid points in the traverse directions are uniformly distributed to match the grid spacing of the HIT domain and avoid interpolation procedures between the two domains. The streamwise grid distribution in the STI domain is determined with the function $x_i = 4 \pi [g(\xi - 0.75) - g(-0.75)]$ , where the normalised grid index $\xi = i/(N_x-1)$ ranges from zero at the first point ( $i=0$ ) to one at the last ( $i = N_x-1$ ). Here, $N_x$ is the number of computational grid points employed along the streamwise direction and $g(x)$ is a stretching function defined as

(4.1) \begin{equation} g(x) = 0.4 (\xi - 0.75) - a_2 [f(\xi - 0.75) - f(- 0.75)], \end{equation}

where $a_2$ is a constant numerically determined to obtain a computational grid of the correct size. In the previous expression, the auxiliary function $f$ is defined as $f(x) = x + \log [\cosh (x b_1)]/b_1$ with $b_1 = (1 + a_2^2)^{3/2}/0.1$ . This mapping concentrates approximately 75 % of the $1600$ streamwise grid points within the first 40 % of the computational domain with almost uniform grid spacing. In all cases, the resulting STI grid spacing yields a minimum ratio of the local Kolmogorov length scale $\ell _k = (\nu ^3/\epsilon )^{1/4}$ to the grid spacing, $\ell _k/\Delta x_i \in [1.31, 1.54]$ , where $\nu$ denotes the kinematic viscosity and $\epsilon$ the TKE dissipation rate. This resolution has been demonstrated to be sufficient to correctly capture the turbulent dynamics throughout the computational domain, including the postshock region. Further numerical verification and grid-resolution metrics are provided in Appendix C. For the STI computational domain, characteristic boundary conditions are imposed at both the upstream and downstream boundaries (Poinsot & Lele Reference Poinsot and Lele1992). In particular, the supersonic inflow condition enforces the velocity, temperature and pressure fluctuations computed in the turbulence generator, which is evolved concurrently to continuously provide primitive-variable profiles to the STI domain. A subsonic outflow boundary condition, which weakly imposes the far-field static pressure of the flow, is utilised at the downstream boundary of the computational domain. This numerical formulation does not consider additional sponge regions as previous studies have demonstrated that Navier–Stokes characteristic boundary conditions are sufficient to weakly enforce pressure without injecting spurious reflections when sufficiently long computational domains are utilised (Williams et al. Reference Williams, Di Renzo, Urzay and Moin2024). The normal shock is nominally positioned at a distance of approximately $0.5 \pi$ downstream of the inflow, with the outflow back pressure weakly imposed to preserve this nominal shock-wave position. The pre-computed back pressure to be weakly imposed at the outflow, which produces zero mean drifting in the shock location, is determined iteratively following the procedure of Larsson & Lele (Reference Larsson and Lele2009). The shock is considered stationary when its mean drift velocity is at least four orders of magnitude smaller than the incoming mean flow velocity.

4.2. Compressible homogeneous isotropic turbulence with linear forcing

The turbulence generator consists of a triply periodic cubic computational domain with dimensions $2 \pi \times 2\pi \times 2\pi$ where forced HIT is computed. Attaining a statistically stationary state requires the inclusion of a forcing term in the momentum and energy conservation equations. Petersen & Livescu (Reference Petersen and Livescu2010) demonstrated that forcing for statistically stationary compressible isotropic turbulence must be done considering the dilatational and solenoidal contributions of the forcing. This parametrisation is crucial to prevent unbounded growth of dilatational to solenoidal kinetic energy $\eta = K_d/K_s$ , and thus to control the associated dissipation ratio $\epsilon _d/\epsilon _s$ .

The HTR solver implements the method outlined in Petersen & Livescu (Reference Petersen and Livescu2010): a linear forcing term that distinguishes between the dilatational and solenoidal components of the velocity fluctuations, $\boldsymbol{v}'' = \boldsymbol{v}_{\boldsymbol{s}}'' + \boldsymbol{v}_{\boldsymbol{d}}''$ , by using the Helmholtz decomposition, namely $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{\boldsymbol{s}}'' = 0$ and $\boldsymbol{\nabla }\times \boldsymbol{v}_{\boldsymbol{d}}'' = 0$ . Two alternative approaches exist to prescribe the target turbulent state: fixing the dissipation ratio $\epsilon _d/\epsilon _s$ or fixing the dilatational-to-solenoidal kinetic energy ratio $\eta = K_d/K_s$ . In this work, we adopt the latter strategy. The forcing appearing on the right-hand side of the momentum equations takes the form

(4.2) \begin{equation} \mathcal{F}_{i} = {C_s}\sqrt {\rho }\,{w}_{s, i} + {C_d}\sqrt {\rho }\,{w}_{d, i}, \end{equation}

where ${{w}}_{s, i}$ and ${{w}}_{d, i}$ are the solenoidal and dilatational components of the density-weighted velocity ${{w}}_{i} = \sqrt {\rho }u_{i}$ (Kida & Orszag Reference Kida and Orszag1990a ). The coefficients $C_s$ and $C_d$ are dynamically determined using a differential control inspired by Bassenne et al. (Reference Bassenne, Urzay, Park and Moin2016) to force the ratio of dilatational-to-solenoidal energy, namely

(4.3) \begin{align} C_s = \frac {\epsilon _{s} - G \big (K_s - K_{s,t}\big )/t_0}{2K_s}\qquad \text{and}\qquad C_d = \frac {\epsilon _{d} - \left\langle {P\frac {\partial {u_k}}{\partial {x_k}}} \right\rangle - G \big (K_d - K_{d,t}\big )/t_0}{2K_d}, \end{align}

where $\epsilon _{s}$ and $\epsilon _{d}$ are the solenoidal and dilatational dissipation rates, respectively, and $K_s$ , $K_d$ represent the corresponding kinetic energies. In (4.3), $G = 500$ is a dimensionless gain parameter that modulates the intensity of the corrective forcing, $t_0$ is the characteristic adaptation time of the forcing that is set equal to the integral turnover time, while $K_{s,t}$ and $K_{d,t}$ are the target solenoidal and dilatational kinetic energies, respectively. The corresponding forcing added to the right-hand side of the total-energy equation is therefore

(4.4) \begin{equation} \mathcal{F}_{\rho E} = \mathcal{F}_{i}{u_i} - \langle \mathcal{F}_{i}{u_i}\rangle \end{equation}

to maintain a statistically stationary mean internal energy.

The numerical investigation is conducted at a constant shock intensity of ${\mathcal M}_1 = 5$ , with variations introduced solely in the upstream turbulence characteristics and the thermochemical properties of the gas. Regarding the former, the study focuses on the influence of upstream turbulence intensity, characterised by the turbulent Mach number $\mathcal{M}_t \approx u_{\textit{rms}} \sqrt {3} / \langle c_1\rangle$ , considering two cases: ${\mathcal M}_t = 0.2$ and ${\mathcal M}_t = 0.4$ . Additionally, the contribution of dilatational motions to the upstream kinetic energy is examined via the parameter $\eta$ , which is varied up to a maximum value of $\eta = 0.1$ . Five different conditions are analysed under the assumption of a calorically perfect gas. These numerical configurations are then repeated for the case of a thermally perfect gas, except the $\eta = 0.001$ case, which corresponds to a nearly solenoidal flow.

Table 1 summarises the full range of HIT conditions and corresponding flow statistics. The spatial resolution, quantified by the average Kolmogorov length scale to grid spacing ratio $\ell _k/\Delta x \in [1.79, 2.87]$ , satisfies standard DNS criteria $\ell _k/\Delta x \geqslant 0.5$ for resolving bulk dissipation statistics (Moin & Mahesh Reference Moin and Mahesh1998). As expected, stronger turbulence intensity and higher compressibility lead to increased velocity and thermodynamic fluctuations, reflected in higher root-mean-square (r.m.s.) values of velocity, density, temperature and pressure. For instance, the normalised pressure fluctuation $p_{\textit{rms}} / (\langle \rho _1\rangle \langle c_1\rangle ^2)$ increases from 0.014 in case C5 to 0.128 in C4. The thermally perfect gas cases (G1–G4) show similar statistics to their calorically perfect counterparts (C1–C4), since the upstream adiabatic index $\gamma$ is almost constant. The parameter $\chi$ , defined in (3.9) as the ratio of entropic-density perturbations and rotational velocity disturbances, remains small across all configurations, indicating that entropic fluctuations contribute negligibly to the upstream turbulence. Consequently, the density and temperature fluctuations approximately satisfy isentropic acoustic scaling. Specifically, the relations $\rho _{\textit{rms}} / \langle \rho _1\rangle \approx p_{\textit{rms}} /( \langle \rho _1\rangle \langle c_1\rangle ^2)$ and $T_{\textit{rms}} / \langle T_1 \rangle \approx (\gamma - 1) \rho _{\textit{rms}} / \langle \rho _1\rangle$ hold across all cases, consistent with the expected behaviour of nearly isentropic acoustic modes in the limit of weak entropy fluctuations (Kovásznay Reference Kovásznay1953).

Table 1. Summary of the numerical simulation set-up for the HIT generator.

Under these conditions, the r.m.s. pressure fluctuation can be approximated by the dilatational velocity component as $p_{\textit{rms}}/ (\langle \rho _1\rangle \langle c_1\rangle ^2) \approx u_{d, {rms}} / \langle c \rangle$ . Expressing $u_{d, {rms}}$ in terms of the turbulent Mach number and the dilatational-to-solenoidal TKE ratio, leads to the scaling

(4.5) \begin{equation} \dfrac {p_{\textit{rms}}}{\langle \rho _1\rangle \langle c_1\rangle ^2} \approx \mathcal{M}_t\sqrt {\dfrac {\eta }{\eta + 1}}. \end{equation}

This expression, reformulated here in terms of $\eta$ , provides a direct relation between normalised pressure fluctuations and the energy content of the dilatational modes in compressible HIT, recovering the same scaling trend reported in previous studies (Sarkar et al. Reference Sarkar, Erlebacher, Hussaini and Kreiss1991; Wang et al. Reference Wang, Gotoh and Watanabe2017b , Reference Wang, Wan, Chen, Xie and Chen2018b ; Donzis & John Reference Donzis and John2020).

While the present results reflect predominantly acoustic fluctuations, other configurations – such as those with prescribed entropy perturbations (Mahesh et al. Reference Mahesh, Lele and Moin1997; Jamme et al. Reference Jamme, Cazalbou, Torres and Chassaing2002) – can exhibit distinct thermodynamic responses across the shock. Additional statistical measures related to vorticity and velocity divergence are discussed subsequently.

As an illustrative example, figure 6 shows the interplay between the turbulent Mach number $\mathcal{M}_t$ and the dilatational-to-solenoidal TKE ratio $\eta$ in shaping the topology of three-dimensional HIT. Since this corresponds to upstream conditions – i.e. a cold gas – the influence of vibrational degrees of freedom, represented through the temperature-dependent specific heat function $c_v(T)$ , is negligible and therefore omitted in this preshock turbulence characterisation. The figure shows instantaneous slices at the $x = \pi$ plane for two values of $\mathcal{M}_t = 0.2$ (top row) and $\mathcal{M}_t = 0.4$ (bottom row), and two TKE ratios $\eta = 0.05$ (left column) and $\eta = 0.1$ (right column), corresponding to cases C1-C4. Panel (a) displays the vorticity magnitude $\varpi = |\boldsymbol{\nabla }\times \boldsymbol{v}|$ , and panel (b) shows the velocity divergence $\vartheta = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}$ . To emphasise the most intense features, the upper limit of the vorticity colour map and both limits of the divergence colour map are set by the 99th percentile of each field, i.e. $\varpi _{99} \equiv {P}_{99}(|\boldsymbol{\nabla }\times \boldsymbol{v}|)$ and $\vartheta _{99} \equiv {P}_{99}(|\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}|)$ . This choice ensures that the visualisations capture the most energetic and intermittent structures in the flow while suppressing background fluctuations. It is worth noting that $\vartheta _{99} \sim {3\vartheta _{\textit{rms}}}$ , aligning with the threshold $\vartheta \lt -3\vartheta _{\textit{rms}}$ commonly used to detect shocklets in compressible turbulence (Samtaney, Pullin & Kosović Reference Samtaney, Pullin and Kosović2001). Although we do not explicitly identify shocklets here, Samtaney et al. (Reference Samtaney, Pullin and Kosović2001) state that a more complete criterion also requires the presence of a density inflection point, ${\nabla} ^2\! \rho = 0$ .

Figure 6. Instantaneous slices from the three-dimensional HIT at $x = \pi$ plane, for different turbulent Mach numbers $\mathcal{M}_t = [0.2, 0.4]$ and dilatational-to-solenoidal TKE ratios $\eta = [0.05, 0.1]$ , using the calorically perfect model (cases C1–C4). Panel (a) displays the vorticity magnitude, $|\boldsymbol{\nabla }\times \boldsymbol {v}|$ , and panel (b) shows the velocity divergence, $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {v}$ . Colour-bar limits are set by the 99th percentile of the corresponding field in case C4, to highlight the growing intensity and spatial intermittency of turbulent structures with increasing $\mathcal{M}_t$ and $\eta$ .

As expected, increasing $\mathcal{M}_t$ from 0.2 to 0.4 intensifies the vorticity field and broadens the range of turbulent scales. Likewise, increasing $\eta$ leads to more prominent dilatational motions and an associated intensification of the vorticity field. The underlying mechanisms responsible for this amplification have been investigated in previous studies. For instance, Kida & Orszag (Reference Kida and Orszag1990b ) showed that, while baroclinic generation can initiate vorticity inside shocklets and the dominant mechanism of enstrophy amplification near shocks is compression, away from shocks, vorticity is primarily enhanced through the classical vortex stretching mechanism. Supporting this, Wang et al. (Reference Wang, Shi, Wang, Xiao, He and Chen2011) reported that shock-induced compression increases tangential vorticity by up to 50 % in their numerical simulations, and that this enhancement persists across a wide range of scales. The combined effect is especially evident when comparing cases C1 and C4, where both $\mathcal{M}_t$ and $\eta$ are increased.

Figure 7. Normalised probability density distributions (p.d.f.s) of vorticity $\varpi / \varpi _{\textit{rms}}$ (a) and velocity divergence $\vartheta / \vartheta _{\textit{rms}}$ (b) for different turbulent Mach numbers $\mathcal{M}_t = [0.2, 0.4]$ , dilatational-to-solenoidal TKE ratio $\eta = [0.05, 0.1]$ and the caloric perfect model (cases C1–C4). Calculations were computed in the HIT domain and averaged over different snapshots.

Although the vorticity field becomes more intense in absolute terms as $\mathcal{M}_t$ and $\eta$ increase, the normalised p.d.f.s shown in figure 7(a) collapse well across all cases, particularly in the bulk region around the mean. Unlike the visualisation in figure 6, which uses 99th-percentile normalisation to highlight spatial extremes, the p.d.f.s in figure 7 are normalised by the r.m.s. to allow consistent statistical comparison of intermittency across different flow conditions. This collapse indicates a comparable level of relative intermittency despite increasing dilatational content. Small deviations appear as $\mathcal{M}_t$ increases (cases C3–C4), where a slightly thinner tail is observed, suggesting a modest reduction in extreme vorticity events relative to the r.m.s.. Nevertheless, this effect is minor, and the overall similarity across cases suggests that the statistical structure of vorticity fluctuations remains largely unaffected by increasing compressibility, consistent with previous findings that solenoidal turbulence statistics are relatively insensitive to compressibility effects (Erlebacher & Sarkar Reference Erlebacher and Sarkar1993).

In contrast, the dilatational component, shown in figure 6(b), is dominated by sharp, sheet-like structures associated with localised regions of intense compression $ (\vartheta /\vartheta _{99} \lt -1 )$ and blob-like regions of intense expansion $ (\vartheta /\vartheta _{99} \gt 1 )$ . These closely spaced compression–expansion regions are consistent with the dynamics of shocklet-like structures, which naturally emerge in high $\mathcal{M}_t$ turbulent flows due to localised steepening of compressive motions, as previously observed in other studies (Pirozzoli & Grasso Reference Pirozzoli and Grasso2004; Larsson & Lele Reference Larsson and Lele2009; Wang et al. Reference Wang, Shi, Wang, Xiao, He and Chen2011, Reference Wang, Gotoh and Watanabe2017a ; Watanabe, Tanaka & Nagata Reference Watanabe, Tanaka and Nagata2021; Sakurai & Ishihara Reference Sakurai and Ishihara2024). The normalised p.d.f.s of velocity divergence in figure 7(b) further support this observation. The distributions $\vartheta /\vartheta _{\textit{rms}}$ show markedly broader negative tails (compression region) for higher $\mathcal{M}_t$ and $\eta$ . This reflects a significant increase in the spatial intermittency of the dilatational field, highlighting the growing dominance of strong localised compression events. These findings are consistent with the results of Wang et al. (Reference Wang, Shi, Wang, Xiao, He and Chen2011, Reference Wang, Shi, Wang, Xiao, He and Chen2012b , Reference Wang, Gotoh and Watanabe2017a , Reference Wang, Wan, Chen, Xie and Chen2018b ), who report that the effects of shocks and compressibility span across a wide range of scales, with their characteristics varying according to the turbulent Mach number and the level of dilatational energy in the flow. Since these shocklet-like features involve steep gradients over very short length scales, it is important to ensure that the simulation is adequately resolved. In this study, the spatial resolution is sufficient to capture such structures, whose characteristic thickness is typically of the order of the Kolmogorov scale (Samtaney et al. Reference Samtaney, Pullin and Kosović2001).