4.1 Stretched-vortex model

The interaction of the resolved flow with the subgrid flows is governed by the subgrid-scale (SGS) stress and heat flux (3.5). The stretched-vortex model (SVM) approximates these interactions by assuming that the SGS flow within each computational cell can be modelled by an ensemble of stretched, spiral vortex flow structures (Misra & Pullin Reference Misra and Pullin1997; Kosović, Pullin & Samtaney Reference Kosović, Pullin and Samtaney2002). The energy spectrum of the subject spiral vortex flow follows a similar cascade to that of isotropic incompressible turbulence, (Lundgren Reference Lundgren1982)

(4.1) $$\begin{eqnarray}\displaystyle E(k)={\mathcal{K}}_{0}\unicode[STIX]{x1D700}^{2/3}k^{-5/3}\text{e}^{(-2k^{2}\overline{\unicode[STIX]{x1D708}}/(3|\tilde{a}|))}, & & \displaystyle\end{eqnarray}$$

where $\tilde{a}=\tilde{\unicode[STIX]{x1D61A}}_{ij}e_{i}^{v}e_{j}^{v}$ is the strain rate, $\tilde{\unicode[STIX]{x1D61A}}_{ij}$ , aligned with the unit vector along the vortex axis, $\text{}\underline{e}^{v}$ and

(4.2) $$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D61A}}_{ij}=\frac{1}{2}\left(\frac{\unicode[STIX]{x2202}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}\tilde{u} _{j}}{\unicode[STIX]{x2202}x_{i}}\right). & & \displaystyle\end{eqnarray}$$

This study assumes the SGS flow consists of a single spiral vortex with axis, $\text{}\underline{e}^{v}$ , aligned with the principal extensional eigenvector of the strain rate tensor. Vorticity-based orientation models allow for the model to produce backscatter (Kosović et al. Reference Kosović, Pullin and Samtaney2002), or anti-dissipative behaviour, which has been suggested to be significant downstream of a shock (Livescu & Li Reference Livescu and Li2017), but previous implementations of the vorticity alignment model have seen relatively poor results in shock–turbulence LES (Bermejo-Moreno et al. Reference Bermejo-Moreno, Larsson and Lele2010). The SGS stress and heat flux for the single, strain-aligned spiral vortex are,

(4.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{ij}=\bar{\unicode[STIX]{x1D70C}}\tilde{k}^{\prime }(\unicode[STIX]{x1D6FF}_{ij}-e_{i}^{v}e_{j}^{v}), & \displaystyle\end{eqnarray}$$

(4.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle q_{i}^{\text{T}}=-c_{p}\overline{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x0394}x}{2}\sqrt{\tilde{k}^{\prime }}(\unicode[STIX]{x1D6FF}_{ij}-e_{i}^{v}e_{j}^{v})\left(\frac{\unicode[STIX]{x2202}\tilde{T}}{\unicode[STIX]{x2202}x_{j}}\right), & \displaystyle\end{eqnarray}$$

$\tilde{k^{\prime }}=\int _{k_{c}}^{\infty }E(k)\,\text{d}k$

$k_{c}=\unicode[STIX]{x03C0}/\unicode[STIX]{x0394}x$

(4.4) $$\begin{eqnarray}\displaystyle \tilde{k^{\prime }}=\frac{1}{2}\overline{\unicode[STIX]{x1D70C}}{\mathcal{K}}_{0}\unicode[STIX]{x1D700}^{2/3}\left(\frac{2\overline{\unicode[STIX]{x1D708}}}{3|\tilde{a}|}\right)^{1/3}\unicode[STIX]{x1D6E4}\left[-1/3,\frac{2\overline{\unicode[STIX]{x1D708}}k_{c}^{2}}{3|\tilde{a}|}\right], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}$ is the incomplete gamma function. The ${\mathcal{K}}_{0}\unicode[STIX]{x1D700}^{2/3}$ prefactor is computed from a matching to the local resolved velocity structure functions, $\langle \overline{{\mathcal{F}}}_{2}(\text{}\underline{x},\text{}\underline{r})\rangle =\langle (\tilde{\text{}\underline{u}}(\text{}\underline{x}+\text{}\underline{r})-\tilde{\text{}\underline{u}}(\text{}\underline{x}))^{2}\rangle$ . The brackets denote an average value, in this case taken over the surface of a sphere of radius $\unicode[STIX]{x1D6E5}$ . This yields (Voelkl, Pullin & Chan Reference Voelkl, Pullin and Chan2000; Hill & Pullin Reference Hill and Pullin2004)

(4.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}_{0}\unicode[STIX]{x1D700}^{2/3}=\frac{\langle \overline{{\mathcal{F}}}_{2}(\text{}\underline{x},\unicode[STIX]{x0394}x)\rangle }{A\unicode[STIX]{x0394}x^{2/3}}, & \displaystyle\end{eqnarray}$$

(4.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle \langle \overline{{\mathcal{F}}}_{2}(\text{}\underline{x}_{0},\unicode[STIX]{x1D6E5})\rangle =\frac{1}{6}\mathop{\sum }_{j=1}^{3}(\unicode[STIX]{x1D6FF}\widetilde{u_{1}^{+}}^{2}+\unicode[STIX]{x1D6FF}\widetilde{u_{2}^{+}}^{2}+\unicode[STIX]{x1D6FF}\widetilde{u_{3}^{+}}^{2}+\unicode[STIX]{x1D6FF}\widetilde{u_{1}^{-}}^{2}+\unicode[STIX]{x1D6FF}\widetilde{u_{2}^{-}}^{2}+\unicode[STIX]{x1D6FF}\widetilde{u_{3}^{-}}^{2})_{j}, & \displaystyle\end{eqnarray}$$

(4.5c ) $$\begin{eqnarray}\displaystyle & \displaystyle A=4\int _{0}^{\unicode[STIX]{x03C0}}s^{-5/3}(1-\sin (s)/s)\,\text{d}s\approx 1.90695, & \displaystyle\end{eqnarray}$$

$(\unicode[STIX]{x1D6FF}\widetilde{u_{i}^{\pm }}^{2})_{j}=(\tilde{u} _{i}(\text{}\underline{x}_{0}\pm \text{}\underline{e}_{j}\unicode[STIX]{x1D6E5})-\tilde{u} _{i}(\text{}\underline{x}_{0}))^{2}$

$\text{}\underline{e}_{j}$

$j$

4.2 Hybrid stretched-vortex model

It was found that the implementation of the stretched-vortex model discussed in § 4.1 developed aliasing errors in the forced upstream turbulence that resulted in an underestimation of the turbulent kinetic energy dissipation rate downstream of the shock. To remove these aliasing errors in the forced turbulence a hyperviscous stress is introduced into the model (Braun, Pullin & Meiron Reference Braun, Pullin and Meiron2018),

(4.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}_{ij}}{\unicode[STIX]{x2202}x_{j}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(\bar{\unicode[STIX]{x1D70C}}\tilde{k}^{\prime }(\unicode[STIX]{x1D6FF}_{ij}-e_{i}^{v}e_{j}^{v}))+\frac{\unicode[STIX]{x0394}x^{5}\bar{c}}{96\unicode[STIX]{x1D6FC}_{max}}\min (\max (\unicode[STIX]{x1D6FC},0),\unicode[STIX]{x1D6FC}_{max})\bar{\unicode[STIX]{x1D70C}}L_{3}[u_{i}], & \displaystyle\end{eqnarray}$$

(4.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle q_{i}^{\text{T}}=-c_{p}\overline{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x0394}x}{2}\sqrt{\tilde{k}^{\prime }}(\unicode[STIX]{x1D6FF}_{ij}-e_{i}^{v}e_{j}^{v})\left(\frac{\unicode[STIX]{x2202}\tilde{T}}{\unicode[STIX]{x2202}x_{j}}\right), & \displaystyle\end{eqnarray}$$

(4.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle L_{n}[\tilde{u} _{i}]=(-1)^{n}\mathop{\sum }_{j=1}^{3}\frac{\unicode[STIX]{x2202}^{2n}\tilde{u} _{i}}{\unicode[STIX]{x2202}x_{j}^{2n}}, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FC}$

(4.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}(\text{}\underline{x})=\frac{\langle \overline{{\mathcal{F}}}_{2}(\text{}\underline{x},\unicode[STIX]{x1D6E5}_{1})\rangle }{\langle \overline{{\mathcal{F}}}_{2}(\text{}\underline{x},\unicode[STIX]{x1D6E5}_{2})\rangle }\left(\frac{\unicode[STIX]{x1D6E5}_{2}}{\unicode[STIX]{x1D6E5}_{1}}\right)^{2/3}-1. & & \displaystyle\end{eqnarray}$$

The scaling of the structure function ${\mathcal{F}}_{2}(\text{}\underline{x},\unicode[STIX]{x1D6E5})\propto (\unicode[STIX]{x1D700}\unicode[STIX]{x1D6E5})^{2/3}$ (Batchelor Reference Batchelor1953) suggests that $\unicode[STIX]{x1D6FC}$ should approach zero in an average sense in a $k^{-5/3}$ inertial range. The constant $\unicode[STIX]{x1D6FC}_{max}$ enforces numerical stability requirements and is taken as $\unicode[STIX]{x1D6FC}_{max}=10$ for all simulations in this study. The separation scales are $\unicode[STIX]{x1D6E5}_{2}=2\unicode[STIX]{x0394}x$ and $\unicode[STIX]{x1D6E5}_{1}=\unicode[STIX]{x0394}x$ , and $\langle \overline{{\mathcal{F}}}_{2}(\text{}\underline{x},\unicode[STIX]{x1D6E5}_{2})\rangle$ is calculated on a restricted mesh with grid spacing $2\unicode[STIX]{x0394}x$ . The hyperviscosity effectively acts as a physically local alternative to spectrally sharp explicit filters, which are typically implemented using global methods such as high-order compact finite differences that are impractical in the localized domain decomposition used in the numerical method discussed in § 5. This model is referred to as the hybrid stretched vorted model (HSVM).

4.3 SGS model near the shock

Subgrid-scale flow in a computational cell containing a shock is significantly different from the type of smooth turbulence that SGS models are typically constructed to represent. The SGS turbulence model becomes excessively dissipative near a shock, and so the SGS interaction model (4.6) is set to zero where the weighted essentially non-oscillatory (WENO) method is active. Previous studies have found WENO alone to be sufficiently dissipative in this region (Bermejo-Moreno et al. Reference Bermejo-Moreno, Larsson and Lele2010).

The application of the SGS model in meshes refined by adaptive mesh refinement (AMR) near to but not directly on the shock remains a problem, because the cutoff wavenumber $k_{c}$ in (4.4) is not clearly defined in regions of AMR. A convection model is introduced that assumes the unresolved kinetic energy is unchanged by the prolongation of the flow from the coarse mesh to the fine mesh in front of the shock, and tracks the motion of the SGS kinetic energy as the flow passes through the narrow band of AMR cells. The model is detailed in appendix A, as the methodology is specific to this SGS model and flow geometry.

6.1 Subgrid statistics

An advantage of the stretched-vortex SGS model is that the subgrid flow used to construct the SGS stresses can be accounted for during post-processing of the LES. Inclusion of the SGS flows alleviates the need to artificially filter DNS or experimental results when validating the LES, and allows the LES to capture effects that may be obscured at the resolved scale. For instance, the SGS statistics are useful when addressing energy amplifications because turbulent fluctuations transferred to small scales by the shock compression may appear to be dissipated if they are mapped to wavenumbers that are not resolved on the computational mesh of the LES.

The SGS Reynolds stresses follow directly from $\unicode[STIX]{x1D70F}_{ij}$ (4.3), and (4.4) allows reconstruction of the SGS kinetic energy in (6.2). Hill & Pullin (Reference Hill and Pullin2004) derived expressions for radial energy spectra of the SGS flows in transverse planes. The vorticity of the subgrid flow is oriented along the spiral vortex axis, and so the SGS enstrophy may be written

(6.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{ij}^{SGS}=\unicode[STIX]{x1D6FA}_{kk}^{SGS}e_{i}^{v}e_{j}^{v}, & \displaystyle\end{eqnarray}$$

(6.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{ii}^{SGS}=2\int _{k_{c}}^{\infty }k^{2}E(k)\,\text{d}k, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FA}_{ij}=\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{j}$

(6.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}_{ij}^{SGS}={\mathcal{K}}_{0}\unicode[STIX]{x1D700}^{2/3}\left(\frac{2\overline{\unicode[STIX]{x1D708}}}{3|\tilde{a}|}\right)^{-2/3}\unicode[STIX]{x1D6E4}\left(\frac{2}{3},\frac{2\overline{\unicode[STIX]{x1D708}}}{3|\tilde{a}|}k_{c}^{2}\right)e_{i}^{v}e_{j}^{v}. & & \displaystyle\end{eqnarray}$$

The streamwise Taylor microscale, $\unicode[STIX]{x1D706}$ , depends on velocity gradients that are not directly resolved in the simulations but, neglecting compressibility effects in the upstream isotropic turbulence, this scale can be related to the dissipation rate by (Pope Reference Pope2000)

(6.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D706}\approx \sqrt{\frac{15\langle \bar{\unicode[STIX]{x1D708}}\rangle }{\unicode[STIX]{x1D700}}}u_{rms}. & & \displaystyle\end{eqnarray}$$

Although this is expected to provide only an approximation to $Re_{\unicode[STIX]{x1D706}}$ , the stretched-vortex model includes a viscous drop off in the subgrid flows and is thus capable of capturing some Reynolds number effects.

8.1 Summary of simulations

Table 1 summarizes the LES simulations discussed in this study. Upstream quantities such as $M_{t}$ and $Re_{\unicode[STIX]{x1D706}}$ contain contributions from the SGS model and thus are approximations. The computational grids are isotropic such that $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=\unicode[STIX]{x0394}z$ , including within regions of AMR. The given meshes are on the coarsest level, such that the simulations with $N_{x}=1024$ , $N_{y}=256$ and an AMR factor of $4\times$ have an equivalent resolution of $4096\times 1024^{2}$ in the direct vicinity of the shock where AMR is active. All simulations are run for a particle passage time across the length of the computational domain, and then statistics are recorded over a subsequent particle passage time or three large eddy turnover times in the upstream turbulence, whichever is longer. This is not sufficient to achieve full temporal convergence of large-scale quantities, and some statistics, such as the ratio of the streamwise to the transverse Reynolds stress, $R_{11}/R_{tr}$ , are slightly perturbed from the value upstream of the shock that would occur with an infinite time average. Given finite computational resources, fine mesh resolutions and a broad parametric study are considered here to be more useful than long simulation times.

Resources are focused on simulations ranging between $M_{s}=1.5$ and $M_{s}=2.2$ because previous DNS have observed interesting behaviour in the post-shock Reynolds stress anisotropy within this regime. DNS has found that, for shocks weaker than $M_{s}\approx 1.5$ , LIA predicts the post-shock Reynolds stress anisotropy reasonably well even at modest $M_{t}$ , but for $M_{s}>1.5$ DNS finds $R_{11}/R_{tr}>1$ in the far-field post-shock flow and this ratio is approximately constant as $M_{s}$ is further increased (Larsson et al. Reference Larsson, Bermejo-Moreno and Lele2013). Conversely, for $M_{s}>1.5$ LIA predicts that an increasingly large fraction of the turbulent kinetic energy is in the transverse Reynolds stress and by $M_{s}=2.2$ predicts $R_{11}/R_{tr}<1$ in the far field of the post-shock flow. The discrepancy between LIA and the referenced DNS grows at larger $M_{s}$ , but it is this transition that is of particular interest.

Table 1. Summary of simulations. The conditions of the upstream turbulence, $M_{t}$ and $Re_{\unicode[STIX]{x1D706}}$ , contain contributions from the model of the subgrid flows and are thus approximations. $L_{x}$ is the length of the domain in the streamwise direction, $x_{s}$ is the shock position and $k_{0}$ is the wavenumber corresponding to the maximum in the energy spectrum of the upstream turbulence. $N_{x}$ and $N_{y}$ are the resolution of the coarsest mesh level in the streamwise and transverse directions, respectively. AMR gives the factor of refinement applied to the mesh near the shock.

8.2 Mesh sensitivity

The two primary mesh parameters in the subject LES are the coarse mesh spacing and the local refinement level at the shock. Global refinement of the coarse mesh increases the fraction of the turbulent kinetic energy that is captured at the resolved scale and is expected to improve the LES results in smooth turbulent regions away from the shock. The local mesh refinement at the shock acts to reduce artificial dissipation from the shock-capturing scheme, reduces spurious behaviour in the entropy field generated by the shock-capturing scheme (Lee et al. Reference Lee, Lele and Moin1997) and allows the shock curvature to be resolved (Lee et al. Reference Lee, Lele and Moin1997; Garnier et al. Reference Garnier, Sagaut and Deville2002).

Figure 2. LES behaviour for $M_{s}=1.5,M_{t}=0.06,Re_{\unicode[STIX]{x1D706}}=500$ on a $96\times 64\times 64$ base mesh with increasing local mesh refinement at the shock. Lines with symbols show the total statistics, including the contribution from the SGS model. (a) Streamwise Reynolds stress, (b) transverse Reynolds stress, (c) density-specific volume correlation $b=-\langle \unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D708}^{\prime }\rangle$ , (d) squared vorticity trace, $\unicode[STIX]{x1D6FA}_{ii}=\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}$ . Dashed line (▵): uniform grid. Dot-dashed line (○): $2\times$ refinement. Dotted line (▫): $4\times$ refinement. Solid line $(+)$ : $8\times$ refinement. Thick solid line: LIA of (7.1).

Test cases are performed to investigate the level of local refinement needed about the shock. The low- $M_{t}$ case considered here is conservative with respect to mesh refinement requirements because as $M_{t}$ is reduced for a given $M_{s}$ , the shock corrugations become smaller and more difficult to resolve. Figure 2 shows the amplification of various statistics of interest through a $M_{s}=1.5$ shock, plotted in the shock normal direction. The statistics are time averaged and spatially averaged in the transverse directions, and contain the contribution from the SGS stretched-vortex flow. The Reynolds stress and vorticity are normalized by linearly interpolating their upstream values to the mean shock position. The density-specific volume correlation $b=-\langle \unicode[STIX]{x1D70C}^{\prime }v^{\prime }\rangle$ , where $v=1/\unicode[STIX]{x1D70C}$ , is normalized by the upstream turbulent Mach number as suggested by scaling in LIA. The density-specific volume correlation $b$ is favoured over $\langle \unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D70C}^{\prime }\rangle$ because of its applications in incompressible mixing of variable-density fluids (Livescu & Ristorcelli Reference Livescu and Ristorcelli2008; Schwarzkopf et al. Reference Schwarzkopf, Livescu, Gore, Rauenzahn and Ristorcelli2011), but for the subject low $M_{t}$ , single component flows $b\approx \langle \unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D70C}^{\prime }\rangle /\langle \unicode[STIX]{x1D70C}\rangle ^{2}$ . The SGS model does not provide an estimate for $b$ at the subgrid scales, and this explains why $b$ tends to be smaller in the LES than predicted by linear analysis. The contribution of the resolved-scale flows to the vorticity is negligible at this mesh resolution and $Re_{\unicode[STIX]{x1D706}}$ . The downstream oscillations in $b$ , and to a lesser extent in the Reynolds stresses, are not an artefact of poor statistical convergence, and previous DNS have seen similar oscillatory behaviour downstream of shocks at low $M_{t}$ (Ryu & Livescu Reference Ryu and Livescu2014). The SGS statistics behave in a non-physical manner when calculated at the shock, resulting in the large spikes near the shock in figure 2, but as discussed in § 4.3 these statistics are not used to calculate SGS terms in the governing equations in the direct vicinity of the shock. These LES are also performed for a $M_{s}=2.2$ shock with similar results.

The streamwise Reynolds stress converges quickly, but the LES requires a factor of $4\times$ AMR at the shock for the rest of the quantities to converge to an acceptable level. Garnier et al. (Reference Garnier, Sagaut and Deville2002) asserts that numerical convergence requires the shock corrugation to be resolved, and this is evaluated in the LES by extracting a triangular mesh of the shock surface from the LES using the shock detection algorithm of Samtaney, Pullin & Kosovic (Reference Samtaney, Pullin and Kosovic2001). The shock does not develop holes at the subject $M_{t}$ and $M_{s}$ , and thus the displacement of the triangularized shock mesh from its mean location can be converted into a continuous, periodic function on a rectangular $64\times 64$ grid. The radial power spectrum of a two-dimensional (2-D) function $\unicode[STIX]{x1D713}(y,z)$ with Fourier transform $\hat{\unicode[STIX]{x1D713}}$ is defined as

(8.1) $$\begin{eqnarray}\displaystyle E_{\unicode[STIX]{x1D713}}^{2D}(k_{r})=\frac{1}{2}k_{r}\int _{0}^{2\unicode[STIX]{x03C0}}|\hat{\unicode[STIX]{x1D713}}(k_{r},\unicode[STIX]{x1D703}_{k})|^{2}\,\text{d}\unicode[STIX]{x1D703}_{k}, & & \displaystyle\end{eqnarray}$$

where $k_{r}$ and $\unicode[STIX]{x1D703}_{k}$ are the radial and azimuthal wavenumbers. The time-averaged radial power spectra of the shock displacement is given in figure 3, plotted against the prediction of LIA applied to (7.1) with $k_{c,LIA}=1024$ . LIA predicts that the shock displacements scale as a function of the energy spectrum as $E(k)/k^{2}$ (Mahesh Reference Mahesh1996), yielding an initially flat spectrum that transitions into a $k^{-11/3}$ slope. The uniform-grid simulation significantly underestimates the curvature of the shock, but the LES exhibits a trend towards $k^{-11/3}$ behaviour as AMR refines the mesh at the shock. The shock displacement and evolution of streamwise statistics both exhibit an acceptable degree of invariance with AMR level for refinement factors of at least $4\times$ .

Figure 3. Time-averaged radial spectrum of shock displacement, for $M_{s}=1.5$ , $Re_{\unicode[STIX]{x1D706}}\approx 500$ and $M_{t}\approx 0.06$ on a $96\times 64\times 64$ base mesh. The symbols show the results from LES with varying levels of AMR at the shock. (▵): uniform grid. (○): $2\times$ refinement. (▫): $4\times$ refinement. $(+)$ : $8\times$ refinement. Dashed line: LIA of (7.1), which follows a $k^{-11/3}$ slope in the inertial range.

The LES approaches the results of LIA as the upstream $M_{t}$ is reduced, as shown in figure 4 for an $M_{s}=2$ shock. The degree of agreement with LIA at $M_{t}=0.03$ is comparable to that observed in shock-captured DNS under the similar flow conditions at $Re_{\unicode[STIX]{x1D706}}\approx 30$ (Tian et al. Reference Tian, Jaberi, Li and Livescu2017), although the ability of the DNS to match LIA was limited by the low $Re_{\unicode[STIX]{x1D706}}$ simulated. Comparison to LIA requires low $M_{t}$ and shock corrugations become smaller and harder to resolve as $M_{t}$ is reduced. This would suggest that additional AMR should be required as $M_{t}$ is reduced to a small value such as $M_{t}=0.03$ , but as shown in figure 4 there is no substantial change in the LES when AMR refinement about the shock is increased in that case. Tian et al. (Reference Tian, Jaberi, Li and Livescu2017) show that shock capturing in DNS is reasonable when the ratio of the Kolmogorov scale to the numerical shock thickness is large, and the corresponding measure in LES would be a ratio of the length scale associated with the smallest resolved upstream eddies to the numerical shock thickness, $\unicode[STIX]{x1D6FF}_{s}$ . Assuming that the smallest upstream eddies are on the scale of the coarse mesh resolution, $\unicode[STIX]{x0394}x_{c}$ , and the numerical shock is a few cells wide on the refined mesh, this will be a fixed ratio that is only dependent on the shock-capturing scheme and the mesh. This ratio is approximately $\unicode[STIX]{x0394}x_{c}/\unicode[STIX]{x1D6FF}_{s}\approx 1.8$ in the LES with $4\times$ AMR, where $\unicode[STIX]{x1D6FF}_{s}$ is taken as the mean velocity difference across the shock divided by the maximum of $|\unicode[STIX]{x2202}\tilde{u} _{1}/\unicode[STIX]{x2202}x|$ , and is similar to the requirement $\unicode[STIX]{x1D702}/\unicode[STIX]{x1D6FF}_{s}>1.4$ found by Tian et al. (Reference Tian, Jaberi, Li and Livescu2017). The tendency of low $M_{t}$ cases to not depend heavily on AMR suggests that the ratio of $\unicode[STIX]{x0394}x_{c}/\unicode[STIX]{x1D6FF}_{s}$ may be a more practical requirement on mesh refinement but, except at extremely low $M_{t}$ , this requirement appears to be met under conditions similar to those required to resolve shock corrugations.

Figure 4. Streamwise (a) and transverse (b) Reynolds stresses averaged over time and cross-sectional planes for a $M_{s}=2.0$ shock and upstream $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $512\times 128\times 128$ base mesh. Solid line: LIA. Dotted line: $M_{t}=0.12$ and $\times 4$ AMR. Dashed line: $M_{t}=0.06$ and $\times 4$ AMR. Dot-dashed line: $M_{t}=0.03$ and $\times 4$ AMR. Pluses: $M_{t}=0.03$ and $\times 8$ AMR.

8.3 Comparison to DNS

DNS is available for validation of the LES at moderate $Re_{\unicode[STIX]{x1D706}}$ . Figure 5 shows the subject LES compared against the $M_{s}=1.5$ shock-captured DNS of Larsson et al. (Reference Larsson, Bermejo-Moreno and Lele2013) for upstream inlet conditions of $M_{t}=0.16$ , $Re_{\unicode[STIX]{x1D706}}=75$ and $k_{0}=6$ . The method used to construct turbulence at the inflow boundary upstream of the shock was different in the DNS and the peak wavenumber is not as well defined as it is in the LES, yielding some difference in the spatial non-dimensionalization of the simulations. The LES of primary interest in this study are performed at higher $Re_{\unicode[STIX]{x1D706}}\geqslant 100$ , but those LES still resolve a comparable fraction of the turbulent kinetic energy relative to this case. The higher resolution LES appears to underestimate kinetic energy dissipation upstream of the shock relative to the DNS when including the contribution from the subgrid continuation, but the behaviour of the resolved scales in the LES at both mesh resolutions are consistent. The subgrid continuation assumes an inertial range, and thus becomes less reliable when the mesh cutoff in the LES is in the viscous scales of the turbulence, particularly at low Reynolds numbers. The agreement of the LES with the DNS is considered reasonable at both mesh resolutions and changes with mesh resolution at the coarsest level are found to be minimal under the resolutions and flow conditions of primary interest in this study, as shown in figure 6.

Figure 5. Streamwise (a) and transverse (b) Reynolds stresses averaged over time and cross-sectional planes for a $M_{s}=1.5$ shock. Inlet conditions are $M_{t}=0.16$ , $Re_{\unicode[STIX]{x1D706}}=75$ and $k_{0}=6$ . The lines show the resolved-scale statistics and the symbols show the total statistics, including the contribution from the SGS model. The LES have a factor of $4\times$ AMR at the shock, and are normalized by their total upstream value interpolated to the mean shock position. Dotted line (▫): LES with a $192\times 128\times 128$ base grid. Dash-dotted line (○): LES with a $96\times 64\times 64$ base grid. Thick solid line: shock-captured DNS of Larsson et al. (Reference Larsson, Bermejo-Moreno and Lele2013) (data from Hickel et al. Reference Hickel, Egerer and Larsson2014).

Figure 6. Streamwise (a) and transverse (b) Reynolds stresses averaged over time and cross-sectional planes for a $M_{s}=1.5$ shock, with upstream conditions of $M_{t}=0.06$ and $Re_{\unicode[STIX]{x1D706}}=500$ . LES are conducted with a factor of $4\times$ AMR at the shock, and results include the contribution from the SGS model. Solid line: LES with $1024\times 256\times 256$ coarse mesh. Circles: LES with $512\times 128\times 128$ coarse mesh. Dashed line: LIA of (7.1).

Figure 7. Streamwise and transverse Reynolds stress amplifications over a shock. LES has upstream flow conditions of $M_{t}\approx 0.06$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ , and is on a $512\times 128\times 128$ coarse mesh with $4\times$ AMR. The DNS was performed by Ryu & Livescu (Reference Ryu and Livescu2014) at $Re_{\unicode[STIX]{x1D706}}=20$ with $M_{t}$ ranging between 0.05 and 0.12 immediately upstream of the shock. Post-shock values are taken at the position of the maximum post-shock streamwise Reynolds stress and include contributions from the SGS model. Streamwise Reynolds stress amplifications are given by the squares (LES), triangles (DNS) and solid line (LIA). Transverse Reynolds stress amplifications are the circles (LES), pluses (DNS) and dashed line (LIA).

A number of additional comparisons of the HSVM-LES to DNS in the canonical shock–turbulence problem and several simplified test problems may also be found in Braun et al. (Reference Braun, Pullin and Meiron2018).

9.1 Shock Mach number effects

The shock’s effect on the turbulent statistics is strongly Mach number dependent. Previous DNS has seen only limited Reynolds number dependency on the Reynolds stress amplification over the shock (Larsson et al. Reference Larsson, Bermejo-Moreno and Lele2013), suggesting that the LES should still agree reasonably well with low Reynolds number DNS regarding this phenomenon. Figure 7 compares the Reynolds stress amplifications from the LES with those seen in $Re_{\unicode[STIX]{x1D706}}=20$ DNS (Ryu & Livescu Reference Ryu and Livescu2014) and LIA. The LIA results show the ratio of the Reynolds stresses far downstream of the shock to their pre-shock state. The post-shock values from the DNS and LES are taken at the position of the maximum streamwise Reynolds stress which occurs in the LES at $k_{0}x\approx 3$ . The upstream turbulent Mach number in the DNS varies between $0.05$ and $0.12$ , and this variation is likely responsible for the lower transverse Reynolds stress amplification in the $M_{s}=2.2$ DNS. Ryu & Livescu (Reference Ryu and Livescu2014) observed that DNS trended towards the results of LIA as $M_{t}$ was reduced, but the lowest $M_{t}$ DNS run still does not identically reproduce the predictions of LIA. The scatter in the Reynolds stress amplifications results from a lack of complete statistical convergence due to the limited number of large eddy turnover times completed over the duration of the LES. The LES agrees well with the streamwise Reynolds stress amplifications predicted by LIA. The transverse Reynolds stress amplifications are lower in LES than in LIA, but the LES still agree reasonably well with the DNS.

The LES dissipation rate and the vorticity are compared to LIA in figure 8. The dissipation rate amplification in LIA is approximated by the predicted amplification in $\langle T\rangle ^{0.76}\langle \unicode[STIX]{x1D6FA}_{ii}\rangle /\langle \unicode[STIX]{x1D70C}\rangle$ , as the LIA analysis does not explicitly consider viscosity. The LES results are sampled immediately downstream of the shock, and small differences in the position where data are taken, as well as minor differences in the application of AMR and WENO between runs, introduce a small amount of scattering in the LES results. LIA predicts that only the transverse components of the vorticity are amplified, whereas the vorticity in the SGS flows in the LES remains nearly isotropic. At the subject Reynolds numbers the downstream vorticity field is expected to return to isotropy very quickly (Larsson et al. Reference Larsson, Bermejo-Moreno and Lele2013), and so this is not considered an important discrepancy.

DNS has observed vorticity amplifications that converge towards LIA more quickly than the Reynolds stresses as $M_{t}$ is reduced (Ryu & Livescu Reference Ryu and Livescu2014). The LES and LIA amplifications in the vorticity and dissipation agree reasonably well for $M_{s}\leqslant 2.2$ , but at higher shock Mach numbers the post-shock dissipation rate appears to be consistently lower than LIA. The amplification in these quantities is expected to be difficult to capture in LES, because they are inherently high-wavenumber phenomena which are almost entirely modelled by the SGS model for the subject flow conditions. The rapid change in length scales across the shock, subject to a fixed cutoff wavenumber in the LES, suggests that LES may have difficulty modelling these small-scale phenomena at large $M_{s}$ . Comparisons between LES and LIA are comparisons of two different models, neither of which represents an exact solution, but subsequent LES are primarily focused on conditions with $M_{s}\leqslant 2.2$ where the LES shows better agreement with LIA.

Figure 8. Squared vorticity trace $\unicode[STIX]{x1D6FA}_{ii}=\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}$ and kinetic energy dissipation rate immediately downstream of the shock, normalized by their upstream values. The LES conditions are the same as in figure 7. The LIA amplification in dissipation rate is approximated by the amplification in $\langle \unicode[STIX]{x1D707}\rangle \langle \unicode[STIX]{x1D6FA}_{ii}\rangle /\langle \unicode[STIX]{x1D70C}\rangle$ . Vorticity amplifications $\langle \unicode[STIX]{x1D6FA}_{ii}\rangle /\langle \unicode[STIX]{x1D6FA}_{ii}\rangle _{0}$ are the circles (LES) and dashed line (LIA). Dissipation rate amplifications are given by the triangles (LES) and solid line (LIA).

9.2 Reynolds number effects

The peak streamwise Reynolds stress downstream of the shock has been observed to decrease when computed in DNS with Reynolds number increasing from $Re_{\unicode[STIX]{x1D706}}=40$ to $Re_{\unicode[STIX]{x1D706}}=70$ , although the change was very small (Larsson et al. Reference Larsson, Bermejo-Moreno and Lele2013). Higher peak Reynolds stresses are observed at higher Reynolds numbers over the range of $Re_{\unicode[STIX]{x1D706}}=10{-}45$ (Ryu & Livescu Reference Ryu and Livescu2014), but, likewise, the magnitude of the Reynolds number effects was limited.

The streamwise Reynolds stress amplifications from the LES at $M_{s}=1.5$ are plotted in figure 9(a) over the range of $Re_{\unicode[STIX]{x1D706}}=20{-}2500$ . The results agree with the DNS in that there is initially a small increase in the amplification of $R_{11}$ with increasing $Re_{\unicode[STIX]{x1D706}}$ over the range $Re_{\unicode[STIX]{x1D706}}=20{-}40$ , and there is a minor decrease between the $Re_{\unicode[STIX]{x1D706}}=40$ and $Re_{\unicode[STIX]{x1D706}}=70$ cases. At $Re_{\unicode[STIX]{x1D706}}>100$ there is no discernible effect of Reynolds number, consistent with the transition to fully developed turbulence observed in many flow geometries near that Reynolds number (Dimotakis Reference Dimotakis2000). The decrease in the amplification of $R_{11}$ over the range of $Re_{\unicode[STIX]{x1D706}}=40{-}100$ , followed by a levelling off with no further Reynolds number effects, agrees qualitatively with the $M_{s}=1.05$ experiments of Kitamura et al. (Reference Kitamura, Nagata, Sakai, Sasoh and Ito2017), although the amplitude of the observed changes was much larger in that study.

The simulations in figure 9(a) are conducted under conditions that range from DNS, at $Re_{\unicode[STIX]{x1D706}}=20$ , to LES with almost no resolved viscosity at $Re_{\unicode[STIX]{x1D706}}=2500$ . Figure 9(b) shows the fraction of the streamwise Reynolds stress that is contained at the subgrid level in the simulations, measured at the location of the peak in the downstream $R_{11}$ . The downstream values are used instead of the upstream conditions because a greater fraction of the kinetic energy is held in the subgrid scales downstream of the shock. The changes in the LES at various $Re_{\unicode[STIX]{x1D706}}$ do not appear to be directly correlated with the fraction of the kinetic energy that is resolved on the computational mesh, suggesting that these trends are not merely a result of numerical resolution.

Figure 9. Amplification in the streamwise Reynolds stress over a $M_{s}=1.5$ shock, and the fraction of that Reynolds stress held in the subgrid scales downstream of the shock. Downstream quantities are sampled at the position of the peak streamwise Reynolds stress. The dashed line has an upstream $M_{t}\approx 0.06$ and the dotted line has $M_{t}\approx 0.18$ . The squares are LES done on a $512\times 128\times 128$ coarse mesh and the circles are done on a $1024\times 256\times 256$ coarse mesh. Both meshes use $4\times$ AMR.

9.3 Density fluctuations

Figure 10 shows radial spectra of density fluctuations and co-spectra of velocity–density fluctuations in the LES. Density fluctuations upstream of the shock are driven by the low turbulent Mach number and are considerably smaller than the fluctuations generated by the shock. There is an increase in the spectra near the mesh cutoff wavenumber resulting from aliasing errors, but both spectra display a $k^{-5/3}$ slope in the inertial range.

Figure 10. Radial power spectra of density fluctuations (a) and co-spectra of density–velocity fluctuations (b) for a $M_{s}=1.5$ shock, with upstream conditions of $M_{t}\approx 0.18$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $1024\times 256\times 256$ base grid with $4\times$ AMR. In (b), the lines show the co-spectrum of $\overline{\unicode[STIX]{x1D70C}}\tilde{u} _{1}$ , and the symbols show the co-spectrum of $\overline{\unicode[STIX]{x1D70C}}\tilde{u} _{2}$ . Solid line (○): $k_{0}x=-1.34$ . Dotted line $(+)$ : $k_{0}x=5.57$ . Dashed line (▵): $k_{0}x=28.3$ .

The downstream density fluctuations arise from a mix of acoustic and entropy modes generated by the shock. LIA predicts that at low $M_{s}$ acoustic modes account for most density fluctuations, but for $M_{s}>1.65$ the entropy modes become dominant (Lee et al. Reference Lee, Lele and Moin1997), at which time the entropy modes grow rapidly with $M_{s}$ . This behaviour may be observed in the LES, as shown in figure 11. The density fluctuations initially remain approximately constant with $M_{s}$ and contain a region of rapid, smooth decay near the shock resulting from the relaxation of exponentially decaying pressure modes. At higher $M_{s}\geqslant 2.2$ , the density amplification begins to grow significantly and the LES reproduces the transition from a smoothly decaying profile to an sharp jump, consistent with the transition to a regime where exponentially decaying pressure modes are comparatively small relative to entropy fluctuations. Velocity dilatation appears in the transport equations for $a_{i}$ and $b$ (Livescu et al. Reference Livescu, Ristorcelli, Gore, Dean, Cabot and Cook2009), and drives behaviour such as the dissipation of $b$ , $\unicode[STIX]{x1D700}_{b}=\overline{v^{\prime }\unicode[STIX]{x2202}u_{k}^{\prime }/\unicode[STIX]{x2202}x_{k}}$ , (Schwarzkopf et al. Reference Schwarzkopf, Livescu, Gore, Rauenzahn and Ristorcelli2011). In incompressible variable-density turbulence these correlations are the result of dilatation produced by molecular mixing of different densities of fluid and, despite weak compressibility effects at the subject $M_{t}$ , this appears to still be the case in the LES. The correlation $\overline{v^{\prime }\unicode[STIX]{x2202}u_{k}^{\prime }/\unicode[STIX]{x2202}x_{k}}$ calculated at scales that are well resolved in the LES is small compared to the observed $\unicode[STIX]{x1D700}_{b}$ , which is the result of subgrid mixing modelled by (4.6).

Figure 11. Averaged density-specific volume correlation $b=-\langle \unicode[STIX]{x1D70C}^{\prime }v^{\prime }\rangle$ in LES with upstream $M_{t}\approx 0.06$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $512\times 128\times 128$ base grid with $4\times$ AMR. Solid line: $M_{s}=1.2$ . Dashed line: $M_{s}=1.5$ . Dotted line: $M_{s}=2$ . Dash-dotted line: $M_{s}=2.2$ . Long-dashed line: $M_{s}=2.6$ . The circle and triangle symbols show LIA of (7.1) at $M_{s}=1.5$ and $M_{s}=2.6$ , respectively.

Figure 12 shows the mass-weighted velocity fluctuations $a_{i}=\langle \unicode[STIX]{x1D70C}^{\prime }u_{i}^{\prime }\rangle /\langle \overline{\unicode[STIX]{x1D70C}}\rangle$ downstream of the shock. This term represents the turbulent mass flux and, along with $b=-\langle \unicode[STIX]{x1D70C}^{\prime }v^{\prime }\rangle$ , is a correlation that tracked in some Reynolds-averaged Navier–Stokes (RANS) models (e.g. Schwarzkopf et al. Reference Schwarzkopf, Livescu, Baltzer, Gore and Ristorcelli2016). Unlike the Reynolds stresses, which have not consistently displayed a return to isotropy within the limited domains possible in numerical studies (Larsson & Lele Reference Larsson and Lele2009; Larsson et al. Reference Larsson, Bermejo-Moreno and Lele2013), $a_{i}$ appears to return to the isotropic state $a_{i}=0$ downstream of the shock, although it does so slowly in the low $M_{t}$ cases. Capturing the behaviour of measures such as $a_{1}$ and $b$ is important in simulations of shock–turbulence interactions, and the production of these statistics by the shock in the LES agrees reasonably well with the predictions of LIA.

Figure 12. Averaged mass-weighted velocity fluctuations $a_{i}=\langle \unicode[STIX]{x1D70C}^{\prime }u_{i}^{\prime }\rangle /\langle \overline{\unicode[STIX]{x1D70C}}\rangle$ . The LES is shown by the symbols and has upstream $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $1024\times 256\times 256$ base grid with $4\times$ AMR. Solid line: $M_{s}=1.5$ LIA of (7.1). Circles: $M_{s}=1.5$ , $M_{t}=0.06$ . Pluses: $M_{s}=1.5$ , $M_{t}=0.18$ . Triangles: $M_{s}=2.2$ , $M_{t}=0.18$ .

9.4 Post-shock anisotropy

DNS have found that the Reynolds stresses do not quickly return to isotropy, whereas the vorticity field rapidly relaxes to an isotropic state (Larsson & Lele Reference Larsson and Lele2009; Larsson et al. Reference Larsson, Bermejo-Moreno and Lele2013). The LES agrees with this, and even at high Reynolds number there is no apparent trend in the Reynolds stresses in the LES that would indicate a return to isotropy. All components of the SGS vorticity in the LES are amplified by the shock, whereas LIA predicts that only the transverse components of vorticity are amplified, but this is considered to be reasonable because of the quick post-shock return to isotropy of vorticity in DNS (Larsson et al. Reference Larsson, Bermejo-Moreno and Lele2013). The SGS Reynolds stresses, which numerically govern the interaction of the SGS and resolved scales of the flow in (4.6), take longer to return to isotropy, and the HSVM-LES is capable of capturing the initial post-shock anisotropy of the Reynolds stresses at the SGS scales relatively well relative to the predictions of LIA (Braun et al. Reference Braun, Pullin and Meiron2018).

9.4.1 Transverse spectra

The averaged Reynolds stresses and vorticity are aggregate quantities that only provide information on the largest and smallest scales in the flow, and say little about the behaviour of the intermediate scales within the inertial range of the turbulence. The radial spectra of the streamwise and transverse velocities upstream and downstream of a $M_{s}=2.2$ shock are shown in figure 13, and these spectra provide insight into how kinetic energy is distributed at all scales in the flow. The modest drop off in the spectra near the coarse mesh cutoff wavenumber is, in part, a physical artefact of using 2-D spectra, as will be subsequently discussed.

Figure 13. Time-averaged radial velocity spectra about a $M_{s}=2.2$ shock, for upstream conditions of $M_{t}\approx 0.18$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $1024\times 256\times 256$ base grid with $4\times$ AMR. The symbols show the resolved spectra in planes located at $k_{0}x=-1.47$ (○), $k_{0}x=5.52(+)$ , $k_{0}x=28.1$ (▵). The averaged spectra of the modelled subgrid flows at each plane are given by the solid, dotted and dashed lines, respectively, and these are extended into the resolved scales to show agreement with the resolved spectra.

Taking ratios of these radial spectra, one can define an anisotropy parameter as a function of wavenumber (Livescu et al. Reference Livescu, Ristorcelli, Gore, Dean, Cabot and Cook2009; Chung & Pullin Reference Chung and Pullin2010),

(9.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D712}^{2D}=\frac{E_{11}^{2D}}{E_{11}^{2D}+E_{22}^{2D}+E_{33}^{2D}}-\frac{1}{3}. & & \displaystyle\end{eqnarray}$$

If the spectra in (9.1) were taken in three dimensions, then clearly $\unicode[STIX]{x1D712}=0$ for all wavenumbers in an isotropic field. Computing $\unicode[STIX]{x1D712}$ with radial spectra, as done here, will generally return $\unicode[STIX]{x1D712}^{2D}\neq 0$ even if the flow is isotropic. A given wavenumber in a one-dimensional spectrum contains contributions from larger magnitude three-dimensional wavenumbers (Pope Reference Pope2000), and the same phenomenon occurs in radial spectra. Taking the definition of the radial power spectrum tensor, with $k_{2}=k_{r}\cos (\unicode[STIX]{x1D703})$ and $k_{3}=k_{r}\sin (\unicode[STIX]{x1D703})$ ,

(9.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D60C}_{ij}^{2D}(k_{r})=\frac{k_{r}}{2}\int _{0}^{\infty }\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D719}_{ij}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}k_{1}, & & \displaystyle\end{eqnarray}$$

and the isotropic velocity spectrum tensor for $k=\sqrt{k_{1}^{2}+k_{2}^{2}+k_{3}^{2}}$ ,

(9.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{ij}(k)=\frac{E(k)}{4\unicode[STIX]{x03C0}k^{2}}\left(\unicode[STIX]{x1D6FF}_{ij}-\frac{k_{i}k_{j}}{k^{2}}\right), & & \displaystyle\end{eqnarray}$$

one obtains expressions for the radial spectra of streamwise and transverse velocities components,

(9.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle E_{11}^{2D}(k_{r})=\frac{k_{r}}{4}\int _{0}^{\infty }\frac{E(k)}{(k_{1}^{2}+k_{r}^{2})}\left(1-\frac{k_{1}^{2}}{k_{1}^{2}+k_{r}^{2}}\right)\,\text{d}k_{1}, & \displaystyle\end{eqnarray}$$

(9.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{E_{22}^{2D}(k_{r})+E_{33}^{2D}(k_{r})}{2}=\frac{k_{r}}{8}\int _{0}^{\infty }\frac{E(k)}{(k_{1}^{2}+k_{r}^{2})}\left(2-\frac{k_{r}^{2}}{k_{1}^{2}+k_{r}^{2}}\right)\,\text{d}k_{1}. & \displaystyle\end{eqnarray}$$

Evaluating (9.4) with $E(k)=k^{-5/3}$ yields a constant value of $\unicode[STIX]{x1D712}^{2D}=1/33$ , illustrating that the 2-D anisotropy parameter will be slightly perturbed from zero in isotropic turbulence. In the case of a finite spectrum, such as (7.1), the effects are more pronounced. Figure 14 shows the anisotropy parameter calculated in the isotropic turbulence in front of the shock. An analytical result for the model isotropic spectrum (7.1) with $k_{c,LIA}=128$ is constructed by performing the same decomposition into vorticity modes used in LIA, and binning the resulting modes by radial wavenumber. The cutoff in the model spectrum results in an abrupt rise in $\unicode[STIX]{x1D712}^{2D}$ at high wavenumbers, which is an artefact of high 3-D wavenumber transverse modes being aliased onto low 2-D wavenumbers in the radial spectra when the modes are sampled in a 2-D plane. The LES results are thus truncated at $k\leqslant k_{max}/2$ to focus on the inertial range, as $\unicode[STIX]{x1D712}^{2D}$ is not considered a useful analysis tool at the largest wavenumbers.

Figure 14. Time-averaged anisotropy parameter in isotropic upstream turbulence upstream of a $M_{s}=1.5$ shock, with $M_{t}\approx 0.18$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $1024\times 256\times 256$ base grid. The solid line shows $\unicode[STIX]{x1D712}^{2D}$ calculated under the assumption of isotropy using the model spectrum (7.1), and the dashed line is the result for an infinite $k^{-5/3}$ spectrum. The circles are the LES, given at a plane located at $k_{0}x=-2$ , just upstream of the mean shock position shock.

Figure 15 plots $\unicode[STIX]{x1D712}^{2D}$ computed in cross sectional planes at various streamwise locations downstream of a $M_{s}=1.5$ and a $M_{s}=2.2$ shock. After some transience close to the shock, the LES shows nearly uniform value in the anisotropy parameter with radial wavenumber over the inertial range. $\unicode[STIX]{x1D712}^{2D}$ appears to return to isotropy over the extent of the inertial range as the flow evolves downstream, with higher radial wavenumbers becoming isotropic more quickly. The exception to this is a small numbers of modes near and below $k_{0}$ which do not return to a profile similar to the isotropic upstream flow, and the $M_{s}=2.2$ case initially sees an increase in $\unicode[STIX]{x1D712}^{2D}$ at these low radial wavenumbers as the flow decays. The very largest scales, $k\approx 1$ , show continuously increasing anisotropy because these motions are on the scale of the domain size and may interact with themselves through the periodic boundary conditions of the domain. This is physically unrealistic but these scales of the flow contain only a small amount of energy. The contribution of 3-D effects to the 2-D spectra limits the direct usefulness of $\unicode[STIX]{x1D712}^{2D}$ in determining when specific scales in the flow relax to isotropy, but this parameter will also be used in subsequent sections to connect the results of this LES to model problems that simplify analysis in three dimensions.

Figure 15. Time-averaged anisotropy parameter downstream of a shock, for upstream conditions of $M_{t}\approx 0.18$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $1024\times 256\times 256$ base grid with $4\times$ AMR. The solid line shows $\unicode[STIX]{x1D712}^{2D}$ calculated from the results of LIA applied to the model spectrum (7.1), and the dashed line is the pre-shock result for (7.1). The dotted line is the result for an infinite $k^{-5/3}$ spectrum. The LES results are given by the symbols, and are calculated in cross-sectional planes downstream of the shock. (a) $M_{s}=1.5$ . The LES spectra are taken at (○): $k_{0}x=2.5$ ; (△): $k_{0}x=7.8$ ; ( $+$ ):  $k_{0}x=24.0$ . (b) $M_{s}=2.2$ . The LES spectra are taken at (○): $k_{0}x=2.3$ ; (△): $k_{0}x=7.7$ ; ( $+$ ):  $k_{0}x=23.9$ .

Figure 16 shows the behaviour of the first few wavenumbers of the radial velocity spectra immediately downstream of the $M_{s}=2.2$ shock. These spectra are taken at planes within the region where exponentially decaying pressure modes contribute significantly to the Reynolds stresses, but the radial spectrum of the streamwise Reynolds stresses predicted by LIA and LES still show remarkable agreement. The transverse velocity spectrum predicted in LIA shows qualitative agreement with the LES within the inertial range, but at low radial wavenumbers the transverse fluctuations see no significant amplification in the LES, and even see damping in the lowest radial wavenumbers as the flow progresses downstream from $k_{0}x=0.7$ to $k_{0}x=2.3$ . It is notable that this damping at the lowest radial wavenumbers is not necessarily a large-scale phenomenon, and would also be consistent with a loss of energy in fine-scale transverse velocity modes with wavenumbers approximately normal to the shock plane that have been aliased onto low radial wavenumbers.

Figure 16. Time-averaged radial power spectra of the streamwise and transverse velocities downstream of a $M_{s}=2.2$ shock, for upstream conditions of $M_{t}\approx 0.18$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $1024\times 256\times 256$ base grid with $4\times$ AMR. The symbols show the LES, and the lines show LIA applied to the model spectrum (7.1). Solid line (○): $k_{0}x=-1.5$ upstream spectra. Dotted line (▵): $k_{0}x=0.7$ . Dashed line (▫): $k_{0}x=2.3$ .

9.4.3 LES with LIA modelling of the shock

It is noted that, while figure 15 shows that the anisotropy parameter over the extent of the inertial range relaxes to a profile similar to that seen in the isotropic upstream turbulence, the relaxation to isotropy appears to be much slower than the development of the Reynolds stress anisotropy, which takes an approximately steady value for downstream positions $k_{0}x>5$ . If the Reynolds stress anisotropy is strongly influenced by the relaxation to isotropy in the inertial range, then the inertial range must relax to isotropy more quickly than suggested by the 2-D anisotropy parameter (9.1).

Describing the scale dependency of anisotropy with respect to 3-D wavenumber is complicated by the shock and inhomogeneous direction in the full shock–turbulence LES, but some estimation of this behaviour may be extracted from homogeneous LES initialized with a post-shock field from LIA. LIA is applied directly to the Fourier transform of a $Re_{\unicode[STIX]{x1D706}}=500$ , $M_{t}=0.18$ , $256^{3}$ forced periodic box LES to produce a 3-D approximation to the far field downstream of the shock (Ryu & Livescu Reference Ryu and Livescu2014; Livescu & Ryu Reference Livescu and Ryu2016), that is anisotropic but remains periodic and homogeneous. This is then used as an initial condition for a decaying $256^{3}$ LES, which acts as a model for the decay of a periodic, anisotropic element of turbulence downstream of a shock. This is referred to as LIA-processed LES, as LIA is used to model the instantaneous interaction of a periodic cube of turbulence with a shock. LES where the shock is explicitly simulated in the computational domain, as done in the other sections of this study, is referred to as shock-processed LES. DNS initialized with the far-field LIA result has been shown to converge towards shock-processed DNS as $M_{t}$ is reduced (Braun et al. Reference Braun, Pullin and Meiron2018), but LIA can alternatively produce a prediction for the near field immediately downstream of a shock. The decay of turbulence initialized with this near field might produce a better agreement to shock-processed DNS, but the near field from LIA is inhomogeneous which complicates boundary conditions and prevents the calculation of 3-D spectra that are the primary motivation for this analysis.

Figure 17. Three-dimensional spectra of streamwise and transverse velocity fluctuations in LIA-processed LES, in which a periodic, anisotropic $256^{3}$ LES is initialized by applying LIA to an isotropic LES with $Re_{\unicode[STIX]{x1D706}}\approx 500$ and $M_{t}\approx 0.18$ . The LIA models a $M_{s}=2.2$ shock. An equivalent downstream position $k_{0}x$ is produced by scaling the simulation time by the mean convection velocity downstream of a stationary shock. Dotted line: $E_{11}(k)$ pre-LIA. Dot-dashed line: $E_{22}(k)$ pre-LIA. Solid line: $E_{11}(k)$ at $k_{0}x=0$ . Dashed line: $E_{22}(k)$ at $k_{0}x=0$ . Circles: $E_{11}(k)$ at $k_{0}x=6.0$ . Triangles: $E_{22}(k)$ at $k_{0}x=6.0$ .

Figure 18. Three-dimensional spectra of streamwise and transverse velocity fluctuations in LIA-processed LES. Pre-LIA flow conditions are the same as in figure 17. Dotted line: $E_{11}(k)$ pre-LIA. Dot-dashed line: $E_{22}(k)$ pre-LIA. Solid line: $E_{11}(k)$ at $k_{0}x=0$ . Dashed line: $E_{22}(k)$ at $k_{0}x=0$ . Circles: $E_{11}(k)$ at $k_{0}x=6.0$ . Triangles: $E_{22}(k)$ at $k_{0}x=6.0$ .

Figure 17 shows the 3-D spectra of the LIA-processed LES before LIA is applied, immediately after the application of LIA, and then after a small amount of post-LIA decay. The time at which these spectra are taken is converted to a spatial location using the mean streamwise velocity downstream of a stationary $M_{s}=2.2$ shock, which allows comparison with the shock-processed LES. The oscillations in the spectra of the transverse velocity fluctuations are a result of the discretization of the upstream field onto integer wavenumbers, as some wavenumbers may not have a significant amount of energy mapped to them by LIA. The transfer of transverse fluctuations towards smaller scales over the shock is clearly visible in the spectra, and the energy content of transverse fluctuations is considerably damped at low wavenumbers. The downstream $k_{0}x$ is quite close to the shock, but the flow has already begun to show a noticeable movement towards isotropy over most of the inertial range. The initial anisotropy in the inertial range varies greatly with shock Mach number, but these scales reliably return to isotropy over the range of $M_{s}$ considered, as shown in figure 18. The flows are initialized using the far-field LIA result, and thus do not contain the exponentially decaying pressure modes that are predicted near the shock. LIA predicts that pressure fluctuations are weighted towards high streamwise wavenumbers downstream of the shock (Ribner Reference Ribner1986a ), and pressure–velocity correlations are often believed to drive the Reynolds stresses towards isotropy. Thus, the actual flow which contains these fluctuations near the shock would be expected to return to isotropy at small scales faster than observed in this test case. The relaxation of these intermediate-scale modes to isotropy may also be affected by the construction of the subgrid model, but it is noted that, unlike many SGS models, the HSVM-LES approach does not assume isotropy in the subgrid flows (Pullin & Saffman Reference Pullin and Saffman1994), which has allowed previous applications of the model to anisotropic flows such as buoyancy-driven turbulence (Chung & Pullin Reference Chung and Pullin2010). The SGS Reynolds stress anisotropy in HSVM-LES of modelled post-shock flows shows qualitatively similar behaviour to that predicted by LIA (Braun et al. Reference Braun, Pullin and Meiron2018), and the shock-processed LES conducted in this study also shows comparable results.

Figure 19 gives the Reynolds stress anisotropy from the same simulation as shown in figure 17. Even though the LIA-processed LES is initialized with the far-field LIA result, which does not contain the exponentially decaying pressure modes that contribute to the rapid rise in $R_{11}$ downstream of a shock, there is still a rapid correction towards $R_{11}>R_{tr}$ within the span of time shown in figure 17. The Reynolds stress anisotropy in the LIA-processed LES agrees well with the shock-processed LES, and both methods level off to an approximately constant Reynolds stress anisotropy beyond $k_{0}x\approx 5$ downstream of the shock. The anisotropy parameter (9.1) calculated in 2-D cross-sectional planes of the LIA-processed LES also exhibits similar behaviour to that seen in transverse planes of the shock-processed LES, as shown in figure 20. Despite difficulties in interpreting the 2-D anisotropy parameter directly, its strong similarity to the model problem of the LIA-processed LES suggests that the scale-dependent relaxation to isotropy in the shock-processed LES proceeds in a similar manner to the LIA-processed LES. The inertial range of the LIA-processed LES relaxes to isotropy quickly, over time scales that would imply that the shock-processed LES has returned to isotropy over much of the extent of the inertial range by the time the Reynolds stress anisotropy has relaxed to its far-field value.

Figure 19. Reynolds stress anisotropy. The symbols show the anisotropy from shock-processed LES with $M_{s}=2.2$ , upstream $M_{t}\approx 0.18$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $1024\times 256\times 256$ base grid with $4\times$ AMR. The line shows LIA-processed LES, in which LES in a periodic $256^{3}$ box is initialized by applying LIA to an LES of homogeneous isotropic turbulence with the same upstream conditions. The temporal decay in the LIA-processed LES is converted to a spatial coordinate using the mean convection velocity downstream of the shock.

Figure 20. Anisotropy parameter calculated from 2-D radial spectra in planes parallel to the shock. The pre-shock flow has $Re_{\unicode[STIX]{x1D706}}=500$ and $M_{t}=0.18$ . The lines show shock-processed LES, in which the shock is explicitly simulated, and the symbols show LIA-processed LES. The temporal decay in the LIA-processed LES is converted to a spatial coordinate using the mean convection velocity downstream of the shock. Solid lines and circles: $k_{0}x\approx 2.3$ . Dashed line and triangles: $k_{0}x\approx 7.7$ . Dotted line: $\unicode[STIX]{x1D712}^{2D}=1/33$ analytical result for an infinite isotropic $k^{-5/3}$ spectrum.

Figure 21. Reynolds stress anisotropy. The symbols show the anisotropy at the location of the peak downstream value of $R_{11}$ in LES with upstream $M_{t}\approx 0.06$ and $Re_{\unicode[STIX]{x1D706}}\approx 500$ on a $512\times 128\times 128$ base grid with $4\times$ AMR. The solid line is the far-field result of LIA. The dashed line is produced by applying LIA to the model spectrum (7.1) and post-processing the results to enforce that all wavenumbers $k/k_{0}>2$ have returned to isotropy downstream of the shock.